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“>  double  emitter  sheath,'  and  3)  surface  emission  ions.  Inclusion 
of  these  'thre^'phenomena  combined  with  'Che.*'  elimination  of 
previous  sheath  approximations  requires  careful  analysis  and 
calculation  of  the  sheath  structure.  It  is  shown  that  the  “Bchrir" 
matching  condition  must  be  generalized  to  insure  that  self- 
consistjtficy  prevails  throughout  the  entire  sheath  and  not  just  at 


the  plasma-sheath  interface.  TuTtheTi?  it  is  ,  shown  that  plasma 
ion  distribution  coming  into  that  sheath  must  have  its  low  energy 
ions  "nut  off"  to  produce  a  self-consistj£it  collisionless  sheath, 
and  that  each  of  these  emitter  sheath  phenomena  reduce  the 
normalized  (  by  plasma  density  )  net  ion  loss  rate  to  the 


emitter . 

Each  of  these  phenomena  also  raise  the  normalized  plasma 
density  adjacent  to  the  emitter.  The  higher  plasma  density  at 
the  emitter  causes  a  greater  increase  in  the  loss  of  hot  plasma 
electron  energy  to  the  emitter  than  the  corresponding  decrease  ir. 

the  loss  of  ionization  energy  (  carried  by  the  ions  )  to  the 

• 

emitter.  Therefore  these  emitter  sheath  phenomena  increase  arc- 
drop.  Within  the  limitations  of  the  current  thermionic  convertor 
formulation,  all  three f  of  these  phenomena  (  which  become 
significant  at  low  currents  )  steepen  the  current -voltage 
characteristic.  1*  At  low  current  densities,  the  present  theory 
shows  that  the  collector  sheath  height  decreases,  resulting  ir.  a 
larger  electron  diffusion  velocity  than  can  be  justified  for  the 
continuum  model  used  in  the  plasma  region.  The  result  of  lower 
performance  at  lower  current  is  in  agreement  with  experimental 
studies,  but  the  limitation  imposed  by  the  formulation  prevent 
theoretical  examination  of  the  lower  current  region  in  which  a 
plateau  of  improved  performance  is  found. 

linclassi^iei 


SECURITY  CLASSIFICATION  OF  THIS  FACE 


'  .  ■.'%  .%  .%  %  AlV.V  %’ 


TABLE  OF  CONTENTS 


DD  Form  1473 . i 

Summary  of  Work . .  .ill 

List  of  Symbols . . . . . v 

1.  Introduction . 1-1 

1.1  The  Cesium  Diode  Thermionic  Convertor . 1-5 

1.2  Formulation  of  the  Problem . 1-10 

1.3  Hie  Emitter  and  Collector  Sheaths  Effects . 1-13 

2.  Thermionic  Convertor  Formulation . II -1 

2.1  General . . . II -2 

2.2  Isothermal . 11-10 

3.  Sheaths . III-l 

3.1  Assumptions . I II -3 

3.2  General  Solution  Condition . 1 1 1  -  7 

3.3  Equations . III-ll 

3.4  Sheath  Solutions . Ill -34 


4.  Isothermal  Solutions . IV-1 

4.1  Effects  of  Ion  Reflection . IV-3 

4.2  Effects  of  Trapped  Ions . IV- 14 

4.3  Effects  of  Surface  Emission . IV-16 

4.4  Comparison  with  Experimental  Work . IV-17 

5.  Non-Isothermal  Solutions . V-l 

5.1  The  Implicit  Computational  Scheme . V-3 

5.2  Non-Isothermal  Results . V-4 


6.  Conclusions  and  Recommendations . 

6.1  Conclusions . 

6.2  Recommendations  for  Further  Work. 

Bibliography 

Appendices 

A.  Asymptotic  Expansions 

B.  Isothermal  Programs 

C.  Non-Isothermal  Programs 


cnrr 


Accession  For 

I>;  7  r  S  OK  MI 
I-VTO  ?A5 

U  T-i.-conced 
-  ttf  ication. 


I.u  s'. r  I  '■‘ivii  i  or./ 

i'  ty  Codes 
C  ii  /.nd/or 

! r  l  s',  i  c  o? 


*"m 

'  .  - 


SUMMARY  OF  WORK 


P 

This  document  contains  a  complete  reproduction  of  the  Ph.D. 
dissertation  of  G.  L.  Main.  The  following  paragraphs  outline  ( 
chapter  by  chapter  )  the  results  of  the  research  carried  out 
under  this  contract. 

|  Chapter  1  develops  the  basic  principles  of  the  thermionic 

convertor  and  its  formulation.  Since  the  formulation  uses  many 
approximations  out  of  computational  necessity,  we  present  an 
}  overview  of  these  approximations  and  the  reasoning  behind  them. 

The  fluid  mechanical  nature  of  the  formulation  yields  some 
fundamental  insight  into  how  the  sheath  affects  the  convertor’s 
performance,  and  gives  some  quantitative  results  on  performance. 

Chapter  2  carries  out  the  formulation  of  the  convertor  for  the 
general  case  and  then  for  the  isothermal  electron  case.  These 
results  are  entirely  conventional.  We  center  on  the  isothermal 
electron  case  for  its  simplicity  in  explaining  emitter  sheath 
effects . 

Chapter  3  on  sheaths  advances  coll is ion less  sheath  theory  to 
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cover  the  sheath  phenomena  under  consideration.  Previous  sheath 
formulations  have  contained  several  mathematical  simplifications 
which  are  incompatible  with  the  study  of  these  sheath  phenomena. 
The  simplifications  are  1)  monoenergetic  ions  from  the  plasma,  2) 
monoenergetic  emitted  electrons,  and  3)  that  the  sheath  potential 
drop  is  large.  Further,  for  these  corrections,  it  is  shown  that 
the  "Bohm"  criterion  for  matching  the  sheath  to  the  plasma  must 
be  generalized. 

Chapter  4  presents  isothermal  results  for  the  thermionic 
convertor  including  the  effects  of  surface  emission,  trapped  ions 
and  reflected  ions.  The  general  result  found  here  is  that, 
contrary  to  intuition,  a  higher  plasma  density  at  the  emitter 
increases  arc-drop.  All  three  emitter  sheath  phenomena  produce  a 
lower  net  ion  loss  rate  and  yield  higher  plasma  density  at  the 
emitter. 

Chapter  5  presents  the  results  of  non-isothermal  calculations 
using  an  implicit  numerical  algorithm  to  confirm  the  isothermal 
results  with  ion  reflection.  However,  detailed  non-isothermal 
calculations  with  trapped  ions  surface  emission  ions  are  not 
carried  out  because  the  combination  of  the  non-isothermal 
algorithm  and  the  full  collisionless  sheath  algorithm  would 
require  large  amounts  of  CPU  time. 

Chapter  6  presents  the  conclusions  of  this  work. 
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CHAPTER  1 


CHAPTER  1:  INTRODUCTION 


1.1  The  Cesium  Diode  Convertor 

1.2  Formulation  for  the  Convertor 

1.3  Emitter  Sheath  Effects 


The  thermionic  energy  convertor  is  perhaps  the  simplest  and 
most  direct  heat  engine  in  existence.  It  converts  heat  directly 
into  electrical  power  by  thermionic  emission.  The  device 
essentially  consists  of  two  electrodes  separated  by  a  gap  ( 
typically  .25  mm  )  containing  cesium  vapor  at  a  low  pressure  ( 

typically  1  torr  ).  One  of  the  two  plates  (  the  emitter  )  is 

heated  by  an  external  source  to  1500  -  1750  K  and  the  other  (  the 
collector  )  is  kept  at  750  -  800  K  by  an  external  sink.  The 

hotter  of  the  two  plates  emits  electrons  thermionically  with  a 
greater  average  velocity  and  a  far  greater  density  than  the 
cooler  plate.  Because  of  the  difference  in  emitted  density  and 
velocity,  a  potential  rise  develops  from  the  emitter  to  the 

collector.  Consequently,  electrical  power  is  generated  directly. 
The  convertor  could  operate  with  a  vacuum  gap;  however  a  vacuum 
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gap  is  impractical  because  an  extremely  small  gap  (  on  the  order 
of  1  micron  or  less  )  would  be  required  to  prevent  the  build  up 
of  a  space  charge  of  electrons  which  limits  current  density.  1 
Cesium  is  introduced  into  the  gap  to  overcome  the  space  charge 
effects  because  cesium  ionizes  easily  (  a  first  ionization  energy 
of  3.d9  eV  ).  The  cesium  diode  convertor,  in  its  ignited  mode,  2 
maintains  a  plasma  electron  temperature  of  approximately  3000  K 
which  supports  the  ionization  of  the  cesium. 


Emitter  sheath  phenomena  are  important  in  thermionic  energy 
convertors  because  the  emitter  sheath  forms  the  boundary 
conditions  for  the  plasma  in  the  gap  and  controls  both  the  ion 
loss  rate  and  the  loss  rate  of  hot  (  3000  K  )  plasma  electrons  to 
the  emitter.  This  thesis  examines  three  expected  emitter  sheath 
phenomena  and  their  effects  on  convertor  performance:  J  1) 
reflection  of  ions  coming  from  the  plasma,  2)  ions  trapped  in  the 
double  emitter  sheath,  and  3)  surface  emission  ions.  Inclusion 
of  these  three  phenomena  combined  with  the  elimination  of 
previous  sheath  approximations  requires  careful  analysis  and 
calculation  of  the  sheath  structure.  It  is  shown  that  the  "bohm" 


1  At  1  micron  approximately  5  amps/cm4  could  pass  through  the 
convertor  under  typical  conditions. 

2  By  ignited  mode,  we  mean  a  plasma  arc  in  which  the  electron 
temperature  is  greater  than  the  temperatures  of  the  electrodes 
bounding  the  plasma. 

3  The  collector  sheath  does  not  have  similar  effects  of  any 
significance  because  its  low  temperature  makes  it  essentially 
non-emitting. 
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matching  condition  must  be  generalized  to  insure  that  self- 
consistancy  prevails  throughout  the  entire  sheath  and  not  just  at 
the  plasma-sheath  interface.  Further,  it  is  shown  that  plasma 
ion  distribution  coming  into  that  sheath  must  have  its  low  energy 
ions  "cut  off"  to  produce  a  self-consistant  collisionless  sheath, 
and  that  each  of  these  emitter  sheath  phenomena  reduce  the 
normalized  (  by  plasma  density  )  net  ion  loss  rate  to  the 
emitter. 

Each  of  these  phenomena  also  raise  the  normalized  plasma 
density  adjacent  to  the  emitter.  The  higher  plasma  density  at 
the  emitter  causes  a  greater  increase  in  the  loss  of  hot  plasma 
electron  energy  to  the  emitter  than  the  corresponding  decrease  in 
the  loss  of  ionization  energy  (  carried  by  the  ions  )  to  the 
emitter.  Therefore  these  emitter  sheath  phenomena  increase  arc- 
drcp.  Within  the  limitations  of  the  current  thermionic  convertor 
formulation,  all  three  of  these  phenomena  (  which  become 
significant  at  low  currents  )  steepen  the  current-voltage 
characteristic.  At  low  current  densities,  the  present  theory 
shows  that  the  collector  sheath  height  decreases,  resulting  in  a 
larger  electron  diffusion  velocity  than  can  be  justified  for  the 
continuum  model  used  in  the  plasma  region.  The  result  of  lower 
performance  at  lower  current  is  in  agreement  with  experimental 
studies,  but  the  limitation  imposed  by  the  formulation  prevent 
theoretical  examination  of  the  lower  current  region  in  which  a 


*  The  limitations  of  the  present  formulation  result  from  the 
asymptotic  division  of  the  plasma  into  a  neutral  plasma  and  a 
collisionless  sheath. 
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plateau  of  improved  performance  is  found. 

The  remainder  of  this  chapter  develops  the  basic  principles  of 
the  thermionic  convertor  and  its  formulation.  Since  the 
formulation  uses  many  approximations  out  of  computational 
necessity,  we  present  an  overview  of  these  approximations  and  the 
reasoning  behind  them.  The  fluid  mechanical  nature  of  the 
formulation  yields  some  fundamental  insight  into  how  the  sheath 
affects  the  convertor's  performance,  and  gives  some  quantitative 
results  on  performance. 

Chapter  2  carries  out  the  formulation  of  the  convertor  for  the 
general  case  and  then  for  the  isothermal  electron  case.  These 
results  are  entirely  conventional.  We  center  on  the  isothermal 
electron  case  for  its  simplicity  in  explaining  emitter  sheath 
effects . 

Chapter  3  on  sheaths  advances  collisionless  sheath  theory  to 
cover  the  sheath  phenomena  under  consideration.  Previous  sheath 
formulations  have  contained  several  mathematical  simplifications 
which  are  incompatible  with  the  study  of  these  sheath  phenomena. 
The  simplifications  are  1)  monoenergetic  ions  from  the  plasma,  2) 
monoenergetic  emitted  electrons,  and  3)  that  the  sheath  potential 
drop  is  large.  Further,  for  these  corrections,  it  is  shown  that 
the  "Bchm"  criterion  for  matching  the  sheath  to  the  plasma  must 
be  generalized. 

Chapter  4  presents  isothermal  results  for  the  thermionic 
convertor  including  the  effects  of  surface  emission,  trapped  ions 
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and  reflected  ions.  The  general  result  found  here  is  that, 
contrary  to  intuition,  a  higher  plasma  density  at  the  emitter 
increases  arc-drop.  All  three  emitter  sheath  phenomena  produce  a 
lover  net  ion  loss  rate  and  yield  higher  plasma  density  at  the 
emitter . 

Chapter  5  presents  the  results  of  non -isothermal  calculations 
using  an  implicit  numerical  algorithm  to  confirm  the  isothermal 
results  with  ion  reflection.  However,  detailed  non-isothermal 
calculations  with  trapped  ions  surface  emission  ions  are  not 
carried  out  because  the  combination  of  the  non-isothermal 
algorithm  and  the  full  collisionless  sheath  algorithm  would 
require  large  amounts  of  CPU  time. 

Chapter  6  presents  the  conclusions  of  this  thesis. 


1 . 1  The  Cesium  Diode  Convertor 

Figure  1.1.1  is  a  schematic  diagram  of  the  cesium  diode 
convertor.  The  emitter  is  heated  externally  to  temperature  T^. 
which  is  typically  1750  K  and  the  collector  is  cooled  to 
temperature  T^,  which  is  typically  750  K.  The  gap  space,  d,  or 
convertor  length,  which  is  typically  .25  mm,  seperates  the 
emitter  from  the  collector.  The  cesium  reservoir,  which  is 
sometimes  imbedded  in  the  collector,  is  kept  at  temperature,  7^, 
to  maintain  the  desired  cesium  pressure  (  typically  1  torr  )  in 
the  gap.  The  electrical  load  is  connected  across  the  emitter  and 
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Figure  1.1.1  The  Cesium  Diode  Convertor 


result  showing  actual  thermionic  convertor  performance  under 
these  conditions.  s  It  should  be  noted  that  at  reasonable  output 

2 

current  densities  (  10  amp/cm  )>  the  output  voltage  is  or.  the 
order  of  0.5  volts.  A  basic  understanding  of  the  thermionic 

convertor  output  is  gained  by  developing  the  ideal  thermionic 

convertor  (  r.c  space  charge  effects  )  in  fig.  1.1.3.  Ir.  the 


*  The  experimental  device  shown  here  is  swept  through  voltage 
from  -0.4  to  0.8  volts  by  an  external  sine  wave  generator  at  60 
Hz.  What  is  shown  as  output  voltage  is  the  applied  voltage.  The 
convertor  is  producing  power  when  the  output  voltage  is  greater 
than  zero. 
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where  is  the  Richardson  emitted  current, 


and  where 

A  =  /2£>  a"pVcm1  k*'.  (4) 

If,  for  instance,  T^  =  1500  K  and  =  2.12  eV  then  =  20 
2 

amp/cm  .  We  have  assumed  the  collector  emits  no  electrons,  and 
therefore  a  lower  collector  work  function  always  improves 

performance.  This  assumption  is  very  good  under  the  usual 

conditions  in  the  thermionic  convertor  (  =  1.6  eV,  =  750  K 

)  because  the  collector  emitted  current  is  10^  times  smaller  then 
the  emitter's  emitted  current. 

The  experimental  results  clearly  show  the  Boltzmann  region, 
except  1)  the  curve  is  steeper  than  exponential,  2)  there  is 
sometimes  a  plateau  at  the  base  of  the  Boltzmann  region,  and  3) 
the  output  voltage  is  lower  than  ideal  because  of  plasma  losses 
known  as  arc-drop.  Previous  work  *  has  shown  that  part  of  the 
steepness  can  be  explained  by  the  ionization  kinetics  of  the 
plasma.  In  this  work  we  show  that  part  of  the  steepness  car.  be 
explained  by  the  sheath  effects.  All  three  of  the  emitter  sheath 
phenomena  considered  here  increase  arc-drop  and  decrease  output 
voltage.  Since  all  three  phenomena  become  more  significant  at 


*  Lawless,  J.L.,  The  Plasma  Dynamics  and  Ionization  Kinetics  of 
The.. ionic  Energy  Conversion,  Ph.D.  Thesis,  Princeton  University, 
1551. 
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low  current  densities,  a  steeper  curve  results.  We  had  also 
wished  to  explain  the  plateau  region  theoretically  because  of  its 
locally  lower  arc-drop.  It  is  believed  that  the  plateau,  which 
sometimes  occurs  at  low  current  density,  results  from  the 
collapse  of  the  plasma  arc  and  is  the  result  of  surface  emission 
as  the  source  of  ions  for  the  plasma.  Unfortunately,  our 
theoretical  calculations  cannot  be  carried  into  this  region 
because  the  formulation  (  the  asymptotic  division  of  the  plasma 
into  a  collisionless  sheath  and  a  neutral  continuum  plasma  region 
)  breaks  down. 


1.2  Formulation  for  the  Convertor 

Best  performance  from  the  thermionic  convertor  is  obtained 
experimentally  when  the  following  empirical  relationship  is 
satisfied: 


p  d  =  10  mil-torr,  (1) 
cs 

where  pcs  is  the  neutral  cesium  pressure  and  d  is  the  gap.  This 
corresponds  to  about  15  ion  mean  free  paths  in  the  gap.  Also, 
under  this  condition,  the  convertor  operates  in  an  ignited  mode 
with  the  plasma  electrons  at  approximately  3000  K  (  the  energy 
necessary  to  maintain  the  plasma  electrons  at  a  higher 
temperature  than  the  emitter  temperature  is  supplied  by  the  arc- 
drop  ) .  At  this  electron  temperature  (  or  lower  temperature  when 
the  convertor  is  net  ignited  such  as  in  the  plateau  ) ,  the 


ionization  is  about  10%. 


This  can  be  seen  from  the  Saha  ecuation 
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for  equilibrium  plasma  density, 


where 


n 

n 


=  neutral  density  (1/cm  ), 


(2) 


i 


=  ion  density  (1/cm  ), 


V, .  =  first  ionization  energy,  and 
f  1 

T@  =  plasma  electron  temperature. 

Under  the  conditions,  p  =1  torr,  T  =  3000  K  and  V 

*cs  ’  e 

eV,  we  have 


fi 


3.69 


nR  =  6.44  X  1015  (1/cm3),  and 
=  8.63  X  1014  (1/cm3). 

Since  the  convertor  is  onny  15  mean  free  paths  long,  the  plasma 
does  not  attain  its  equilibrium  plasma  density  and  recombination 
is  usually  negligible.  Therefore  the  ions  are  generated  by  the 
high  electron  temperature  in  the  gap  and  are  lost  to  the  emitter 
and  collector  by  diffusion.  Surface  emission  of  ions  is 
generally  not  a  factor  in  supporting  the  plasma  as  can  be  seen 
from  the  Saha -Langmuir  equation, 


where 


(3) 


J  .  =  surface  emission  ion  current, 
cs+ 

n  =  neutral  cesium  number  density,  and 
cs 

M  -  cesium  atom  mass, 
cs 

Under  the  previously  used  conditions, 

J  ,=  1.6  X  10  3  amps/cm^. 
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In  the  convertor,  Debye  length. 


which  is  length  scale  for  charge  seperation  in  the  plasma  and 
therefore  the  sheath  length  scaling,  is  small  compared  to  both 
the  electron  and  ion  mean  free  paths, 

A0  K  <<  d  .  (5) 

Therefore,  we  divide  the  convertor  into  a  neutral  plasma  region 
terminated  by  collisionless  sheaths  at  the  electrodes.  Figure 
1.2.1  is  potential  distribution  in  the  convertor  (  definitions  of 
terms  are  detailed  in  chap.  2  ).  Based  on  the  asymptotic 
matching  of  the  neutral  plasma  region  to  the  collisicnless 
sheaths,  the  sheath  results  produce  boundary  conditions  for  the 
neutral  plasma.  The  plasma  region  is  treated  as  fluid  (  with  a 
source  term  for  ionization  )  with  conservation  of  mass,  momentum, 
and  energy.  Conservation  of  mass  requires  a  boundary  condition 
for  the  ion  loss  rate  which  is  found  from  the  sheath  net  ion  flux 
rate.  Conservation  of  momentum  results  in  a  net  change  in 
potential  through  the  plasma  region  which  is  added  to  the  change 
through  the  sheaths.  And  conservation  of  energy  (  which  is 
dominated  by  electron  energy  )  requires  boundary  conditions  or. 
plasma  electron  temperature.  Examination  of  fig.  1.2.1  shows 

that  ions  leaving  the  plasma  for  the  emitter  (  at  0+  )  C2n  be 
reflected  by  the  back  of  the  sheath  if  is  greater  than  zero. 

Ions  entering  the  sheath  are  accelerated  by  the  front  sheath  and 

decelerated  by  the  back  sheath.  Those  ions  entering  the  sheith 
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Figure  1.2.1  Potential  in  the  Convertor 


with  energies  less  than  Axg  reflected  back  through  the  sheath 
into  the  plasma  region  again.  Figure  1.2.2  is  plasma  density  in 
the  convertor  gap  normalized  by  net  current,  J,  and  electron 
speed  of  sound,  c,  for  different  ion  reflection  conditions  at  the 
emitter.  The  higher  curves  are  of  littie  or  no  ion  reflection 
and  the  lower  curves  result  from  ion  reflection. 


1.3  Emitter  Sheath  Effects 
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through  the  plasma  and  therefore  reducing  arc-drop.  However, 
this  is  not  the  case.  While  the  plasma  density  at  the  emitter 
increases  slightly,  plasma  density  at  the  collector  decreases. 
Consequently  total  resistance  increases. 


( 
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CHAPTER  2:  THERMIONIC  CONVERTOR  FORMULATION 


2.1  General  Formulation 

2.2  Isothermal  Formulation 


In  this  chapter  we  develop  the  formulation  for  the  quasi¬ 
neutral  plasma  region  of  the  thermionic  convertor  and  the 
boundary  conditions  for  it.  The  boundary  conditions  contain 
fluxes  of  electrons  and  ions  that  must  be  determined  from  the 
sheaths  and  therefore  the  sheath  is  critical  to  the  formulation. 
Sheaths  are  considered  in  detail  in  the  next  chapter.  The  first 
section  is  the  general  formulation  for  the  non-isothermal  plasma 
electron  case  (  finite  electron  thermal  conductivity  ),  and  the 
second  section  is  specialized  to  the  isothermal  case.  The 
general  formulation  follows  closely  the  notation  of  Lawless  1  and 
the  isothermal  formulation  follows  closely  the  notation  of  Lam. 
*  However,  the  isothermal  formulation  has  been  generalized  to 

1  Lawless,  J.L.,  The  Plasma  Dynamics  and  Ionization  Kinetics  of 
Thermionic  Energy  Conversion,  Ph.D.  Thesis,  Princeton  University, 
1981,  chap.  2. 
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eliminate  the  assumption  of  high  sheaths  which  has  previously 
been  used  to  simplify  the  electron  dynamics.  Lawless  developed 
an  explicit  , unsteady  numerical  scheme  for  the  general  non- 
isothermal  case.  However,  the  explicit  scheme  is  slow  because 
stability  requires  small  time  steps.  In  appendix  C  we  develop  an 
implicit  numerical  scheme  which  eliminates  the  stability 
constraint  and  speeds  up  the  Lawless  scheme  by  the  number  of 
space  grids  squared.  Even  the  implicit  scheme  is  slow  and  it  is 
also  difficult  to  gain  physical  insight  from  direct  numerical 
results.  Therefore  we  develop  the  isothermal  formulation  which 
is  faster  and  easier  to  analyse.  Most  of  the  results  in  this 
thesis  are  based  on  the  isothermal  formulation.  Appendix  B 
contains  the  isothermal  programs . 


2.1  General  Formulation 

In  this  section,  which  follows  the  notation  of  Lawless  ( 
chap. 2  ),  we  develop  the  general  non-isothermal  formulation  for 
the  quasi-neutral  region  of  the  thermionic  convertor.  We  assume 
that  the  convertor  is  one-dimensional  since  convertors  typically 
have  a  gap  of  .25  mm  while  the  plates  are  10  cm  in  diameter. 
Following  Lawless,  we  present  the  conservation  equations,  then 
reduce  them  to  convenient  and  useful  forms.  We  next  present  and 
discuss  the  approximations  to  be  used,  and  resulting  parabolic 
P.D.E.s  (  first  order  in  time  and  second  order  in  space  )  which 

1  Lam,  S.H.,  Preliminary  Report  on  Plasma  Arc-Drop  in  Thermionic 
in  Thermionic  Energy  Converters.  Princeton  Unversity,  1976. 
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are  used  to  carry  out  the  numerical  calculations.  Finally,  we 
discuss  the  boundary  conditions  for  these  equations  whose  flux 
rates  (  electrons  and  ions  )  are  derived  from  the  sheaths. 

The  conservation  equations  we  develop  are:  mass  of  ions  and 
electrons,  momentum  of  ions  and  electrons,  and  energy  of 
electrons . 

Conservation  of  mass  of  ions  and  electrons  is: 

sr.  £a  an  (1) 

at  lx  >  J  T* 

where 

T  =  n  u  =  electron  flux, 
e  e  e 

Ti  =  n^u^  =  ion  flux, 

ue  =  mean  electron  velocity, 

u^  =  mean  ion  velocity, 

nft  -  electron  number  density, 

n^  -  ion  number  density, 

S(n)  =  net  rate  of  ionization. 


Conservation  of  momentum  for  ions  and  electrons  is: 
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where 

m  =  electron  mass, 

M  =  ion  mass, 

Pe  =  electron  pressure, 
p^^  =  ion  pressure, 
q  =  electron  charge, 

*  potential, 

S  ^  =  exchange  of  momemtum  from  neutrals  to 

electrons, 

S  -  exchange  of  momemtum  from  ions  to 

electrons , 

S.a(p)  =  exchange  of  momemtum  from  neatrals  to 

ion. 

Note:  in  this  thesis,  potentials,  $,$,x,  are  defined  such  that 
increasing  potential  repels  electrons.  This  is  in  contrast  to 
the  usual  convention  and  is  done  so  because  electrons  are  of 

principal  interest. 

Conservation  of  electron  energy  is: 


i 


» 


L 


m  j.vm 


* 

r 


(3) 
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2 

e  =  mu  (2  +  3kT  / 2  =  stagnation  energy, 
e  e 

h  =  e  +  Pe/nfi  =  stagnation  enthalpy, 

S  '  =  collisional  and  ionization  energy  source 
term, 

q  =  electron  heat  conduction  flux. 


Other  useful  conservation  equations  can  be  derived  from 
eqs.1,2  and  3.  Also,  from  this  point  forward  we  invoke  the 
quasi-neutral  approximation,  n  =  ng  =  From  conservation  of 

mass  (  eq.l  )  we  have 

=  0  (4) 

where 


31 

9* 


J  =  qCfg  -  =  net  current  density. 


Equation  4  results  from  subraction  of  eqs.l  and  using  the  quasi¬ 
neutral  approximation  to  equate  the  time  derivatives  of  density. 
From  conservation  of  momentum  for  ions  and  electrons,  we  have 


Mc  /  Mn  •Sif' 

kt  D,i 

where 


(5) 


P  =  Pe  +  P^  =  total  plasma  pressure. 

We  have  added  the  two  momentum  equations  and  consequently  have 

eliminated  the  terms,  S  and  n3(q<|>)/3x,  which  represent 

6 1 


( 
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exchange  of  momentum  between  electrons  and  ions  and  the  momentum 
exchanged  with  the  field.  From  electron  energy  conservation,  we 


derive 


&n(e+tV]+&[r.o,+zthle]  = 


+  s 1,1 


by  differentiation  of  a  product.  Finally,  the  electric  potential 
term  can  be  eliminated  from  the  energy  equation  by  subtracting 
the  electron  momentum  equation  multiplied  by  u^: 

+  -  >u(£sf) 


These  three  conservation  equations,  which  are  just  the  first 
three  moments  of  the  Boltzmann  equation,  cannot  be  solved  without 
constitutive  equations  for  pressure,  heat  flux  and  the 
collisional  terms.  In  order  to  produce  the  necessary 
constitutive  relations,  we  assume  that  the  electron  and  ion 
distributions  in  the  plasma  are  near-Maxwellian  and  therefore  we 
can  approximate  (  ideal  gas  laws  ) , 

fa  z  n  fclg, 

f;  =  n^T;  (8) 


where 


.•>  - 


|  I  fc  ^  ^  ^  *  —  0  _  •  ^  I  _  .*  * 
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=  plasma  electron  temperature, 

=  Tn  =  ion  and  neutral  temperature, 

k  =  Boltzmann's  constant. 

Also,  because  of  the  near-Maxwe Ilian  assumption,  we  can  make  the 
first  order  approximation, 

„  _  ,  rr, 

1*  9*  (5> 

where 


kg  =  electron  thermal  conductivity. 

The  previous  two  models  (  electron  heat  flux  and  electron 
pressure  )  are  not  totally  consistent  since  eq.8  assumes  an 
isotropic  Maxwellian  electron  distribution  while  the  existence  of 
an  electron  heat  flux  requires  a  non- isotropic  distribution. 


For  the  collision  terms,  we  assume  momentum  transfer  is 
proportional  to  the  difference  of  mean  species  velocity. 


SeT  = 

Scf  =  -vKifc-n) 

-Mcr, 


where 


v  =  electron-atom  collision  frequency, 
vgi  *  electron-ion  collision  frequency, 
v.  =  ion-atom  collision  frequency. 
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The  quantities  v  ,  v  , ,  v,  and  k  must  be  determined.  Tbis  is 
covered  extensively  by  Lawless.  One  further  assumption  is  used, 
namely,  that  the  substantial  time  derivatives  in  the  momentum 
equations  (  Du/Dt  )  are  small.  This  is  justified  for  two 
reasons:  first,  3u/3t  is  small  because  the  mean  time  between 
collisions  is  on  the  order  of  10  nsec  while  all  characteric  time 
scales  in  the  equations  are  much  greater,  and  second,  u3u/3x  is 
small  when  the  Mach  number  is  small  (  electron  or  ion  ).  The 
first  part  is  always  easily  satisfied.  The  second  part  is  not 
well  satisfied  near  the  electrodes  for  the  ions,  and  not  well 
satisfied  near  the  electrodes  for  the  electrons  when  the  sheath 
heights  are  low. 


We  now  introduce  the  definition  of  ambipolar  flux  and 
resulting  equations: 


where 


>  rrjyma  * 


(ID 


(12) 


The  quantity  T  is  called  ambipolar  flux  and  y.,  and  y  are  called 
mobilities.  With  these  definitions  and  assumptions  we  can  write 
eq.5  (  ion  and  electron  momentum  )  as 


Similarly,  we  can  write  conservation  of  mass  as 


(13) 


9n 


(14) 
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and  conservation  of  energy  as 


1  a  A  (k%)  =  £  /fes )  _  sr  sarj 

z  "  z'e  -57- 

'  njAu  e  '***1 

Two  final  approximations  are  used  to  close  these  equations: 

Sco=  -E.S" 


5^-  fart-  PneZ)  fie. 


The  first  states  that  ionization  is  the  only  volume  source  term 
for  electron  energy  where  is  first  ionization  energy.  The 
second  is  a  model  for  ionization  and  recombination.  The  quantity 
a  called  the  ionization  coefficient,  ft  is  the  recombination 
coefficient  and  NQ  is  the  neutral  density.  Again,  the  choice  of 
this  model  and  its  coefficients  is  covered  in  Lawless. 


Now  that  the  conservation  equations  and  the  constitutive 
relations  are  established,  we  can  derive  boundary  conditions  for 
mass  and  energy.  These  boundary  conditions  contain  fluxes  of 
electrons  and  ions  which  are  obtained  from  the  sheath  results. 
The  boundary  conditions  for  mass  come  from  eq.13, 

r'=  ~*$;l  * 

(18) 


The  quantities  u,  and  J  are  the  fluxes  through  the  sheaths.  The 
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subscripts  0  and  1  denote  the  emitter  and  collector  plasma-sheath 
interfaces  respectively.  The  boundary  conditions  for  electron 
energy  are 


where 

Tg  =  flux  of  electrons  from  the  emitter  into  the 
plasma, 

Tg  —  emitter  temperature, 

Vj.  =  emitter  sheath  height, 

Vc  -  collector  sheath  height. 

We  have  assumed  the  collector  emits  nothing. 

2.2  Isothermal  Formulation 

In  this  section  the  corrected  isothermal  model,  based  on  the 
isothermal  model  of  Lam,  and  including  ion  reflection  at  the 
emitter,  is  developed.  Since  we  encounter  both  low  emitter  and 
low  collector  sheath  heights  as  a  consequence  of  ion  reflection 
and  trapped  ions,  the  assumption  of  Boltzmann  plasma  electron 
distributions  at  the  plasma-sheath  interface  must  be  abandoned. 
At  both  the  emitter  and  collector  the  low  sheaths  return  few 
plasma  electrons,  leaving  the  distributions  largely  one  sided. 
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Furthermore,  at  the  emitter  sheath  emitted  electrons  must  be 
taken  into  account.  Thus  the  ratio  of  electrons  moving  toward 
the  sheath  to  the  total  density  of  electrons  at  the  sheath  edge 
is  not  1/2,  as  in  the  Boltzmann  assumption. 
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where 

X  *  nondimens ional  potential, 
i  ~  potential, 
q  a  electron  charge, 
k  =  Boltzmann's  constant, 

Tg  =  emitter  temperature. 

We  also  use  the  following  terminology  for  various  potentials  in 
the  convertor: 

=  emitter  work  function, 

Ax  =  back  sheath  height, 

AXg  =  reflective  potential, 

Xj.  =  emitter  sheath  height, 

AXp  =  plasma  potential  drop, 

=  arc-drop, 

Xj,  =  collector  sheath  height, 

#£  =  collector  work  function, 

VQut  =  convertor  output  voltage. 

Inspection  of  fig.  2.2.1  yields  immediately  the  following 
relations  which  will  be  useful  later: 
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Vj  =  Ymi  ~  ($t~  ~  iX)  (2) 

Vi  =  (Xc-Xs)  -  .  (3) 

Some  further  preliminary  definitions  are  also  needed.  We  have 
the  Richardson  current  density  of  electrons  from  the  emitter, 

JR  ( =  120  '*P  (-is)  (4, 

The  emitted  current  density  which  crosses  the  emitter  sheath 
potential  peak  into  the  convertor  plasma  region  is 

-r  - 

(5) 

J£  =  JR  ;  AX  ^  O  , 

We  also  define  the  net  current  density  through  the  convertor,  J, 
and  the  normalized  current  density, 

j  =  3/Jg  <6> 

Electron  temperature  is  nondimensionalized  as 

r  =  Te  /rE  (7) 

where  Te  is  the  plasma  electron  temperature  which,  in  this 
section,  is  constant  by  the  isothermal  assumption.  Finally  we 
have  the  thermal  speeds, 

A-  =  xJjtf  (8) 
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(9) 


The  isothermal  formulation  is  developed  from  here  in  the  same 
way  as  the  general  formulation  except  that  we  take  full  advantage 
of  the  isothermal  assumption  by  looking  only  at  the  global 
conservation  equations  instead  of  the  local  ones  used  in  the 
general  formulation.  We  then  assume  that  the  transport 
properties,  collision  frequencies,  and  the  ionization  source 
coefficient  are  constant  across  the  convertor  because  of  the 
isothermal  assumption.  Also  we  find  only  the  steady  state 
solution.  We  carry  out  this  development  by  deriving  the  global 
conservation  equations  for  the  isothermal  case  (  current, 
momentum,  and  electron  energy  )  and  then  reducing  these  to  a  set 
of  three  simultaneous  equations  in  the  variables  t,  x£,  and  x^. 
3  These  equations  are  nonlinear  and  solved  numerically  using  a 
positive  definite  Newton's  method  explained  in  appendix  B. 


First  we  consider  conservation  of  current.  The  collector  is 
assumed  to  emit  nothing,  therefore  at  the  collector  plasma-sheath 
interface  we  have 


CT  = 


0e«,nO) 

z. 


(10) 


where  is  the  fraction  of  the  total  plasma  density  at  the 


3  In  some  cases  the  the  actual  calculations  are  carried  out  using 
different  variables  when  Xg  or  Xj,  are  small  or  zero.  In  the  case, 
for  instance,  of  a  single  ion  repelling  emitter  sheath  we  use  j 
because  Xg  is  zero. 
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collector  sheath  which  is  moving  toward  the  collector  and  n(l)  is 
the  total  plasma  density  at  the  plasma-collector  sheath 
interface.  Because  we  continue  to  assume  that  the  plasma 
electron  distribution  there  is  Maxwellian,  we  can  write  as 


I 


which  takes  into  account  the  plasma  electrons  reflected  by  the 
collector  sheath.  We  still  assume  that  the  plasma  electron 
distribution  coming  into  the  collector  sheath  is  Maxwellian  and 
that  it  does  not  have  any  velocity  shift  because  the  sheath  is 
expected  to  be  electron  repelling.  In  the  limit  of  a  high 
collector  sheath,  =  1/2  and  we  have  a  fully  Boltzmann 

distribution  of  electrons  at  the  collector  sheath  edge.  The 
situation  at  the  emitter  is  more  complex  because  the  emitted 
electrons  must  be  taken  into  account.  We  have  the  back  scattered 
current  density,  JgS»  which  is  the  plasma  electron  current 
density  moving  into  the  emitter, 


T  _  n(0)  <7e  tip  “  X*r  (12) 

w'bs  2.  c 

where  n(0)  is  the  total  plasma  density  at  the  emitter  sheath- 
plasma  interface  and  is  the  fraction  of  total  plasma  density 
at  the  interface  moving  toward  the  emitter. 


Continuity  of  electron  current  demands 

Te=T6S  +  J 


(13) 


which  can  be  written  as 
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7(1+ 

This  can  be  rewritten  using  eqs.3  and  6  as 


The  quantity  can  be  written  as 


where 


VT 


(14) 


(15) 


(16) 


is  the  electron  Mach  number  at  the  emitter.  This  in  just  an 
application  of  eq.13. 


Electron  energy  conservation  is  developed  by  considering 
energy  exchange  with  the  emitter  and  collector  and  energy  lost  to 
ionization.  Power  carried  into  the  plasma  by  emitted  electrons 
is 


(17) 


Power  returned  to  the  emitter  is 

r,s=  (j'-VCs-T+b+u)*?  ■ 

Power  flowing  into  the  collector  is 

Pt=  7 ( 2r  +  &  +  t.x+Vi) &  . 


(18) 


(19) 


Ionization  power  loss  is 
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fie n  Vf ^£' 


(20) 


where  J,  is  the  total  ion  current  into  both  the  emitter  and 
ion 

collector,  and  is  the  first  ionization  energy.  Conservation 


of  electron  energy  is 

P,= 


or\ 


(21) 


this  can  be  reduced  to 

i~ijVj  -±j)v4f  (22) 

where  =  ^ion^E"  t*ie  ignited  mode  t  is  generally  about  2  ( 
Tg  =  1500  K  and  Tg  =  3000  K  ),  consequently  the  arc-drop,  V^,  is 
negative.  In  other  words  the  high  plasma  electron  temperature  is 
generated  by  resistance  heating. 


Finally,  we  consider  electron  and  ion  momentum.  From  electron 
momentum  conservation,  we  find  the  potential  drop  in  the  plasma 
region.  By  adding  the  electron  and  ion  momentum  equations  as  in 
the  general  case,  we  find  our  diffusion  equation  and  boundary 
conditions  to  which  the  sheaths  contribute  flux  terms.  When  we 
introduce  the  ionization  source  term  into  this,  we  have  the 
complete  formulation.  Electron  momentum  conservation  is 

dlk  _  Qm.  r»r)Ue 


(23) 


where  X  is  electron  mean  free  path.  Using  p  ■  nkT  and  J  = 
©  ©  © 

qnu^,  we  can  rearrange  eq.23  into 


7  " 
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This  can  be  further  reduced  by  dividing  by  Jg  and  using  £  =  x/d 
where  d  is  the  convertor  gap  thickness: 

i  -  -X  A*  /  frdn^„dl\ 

J  *  ttffA  ft'nft)-  (25) 

Integration  of  this  equation  from  the  emitter  sheath  interface  to 
the  collector  sheath  interface  yields 


where 


*■  %  ^  ( &■»  • 


The  quanity  R  is  the  normalized  plasma  resistance. 


The  ion  and  electron  momentum  equations  can  be  written 

kL&l  =  -  m'lU'Q'. 

dK  &  dx 

jcT-ia.=z  init-rfnu.Q,- 

d  X  v  dx  ^ . 

where  is  ion  mean  free  path  and  a^  is  ion  thermal  speed, 


A.  =  I /ML 

/  7Ttf  * 


Addition  of  eqs.28  yields 

which  is  ambipolar  diffusion.  Equation  29  is  differentiated  to 


become 


<30) 
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At  this  point  we  assume  recombination  is  negligible  and  the 
ionization  source  term  is 


(31) 


Using  eqs.31  in  eq.30  yields 


djr\ 


t 


(32) 


Equation  29  taken  at  the  boundaries  of  the  plasma  (  at  the 
emitter  and  collector  sheath  interfaces  )  forms  the  plasma 
boundary  conditions 


(33) 


where 


_  -  d 


h  MTt+£re  (  Ac  1 


J 


Equation  32  is  written  as 


din  +  A'Mn-o 

dS1- 


where 


Kcr)  = 


(34) 


(35) 


(36) 


where  A(t)  is  the  ionization  coefficient  and  is  found  from 
consideration  of  ionization  kinetics.  Its  solution  for  n  is 


n(S)  =  Bsin(A?+  O 


(37) 
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where  B  and  C  are  constants  of  integration  and  A=A(t).  The 
quantities  and  p^,  which  are  the  boundary  conditions  for 
eq.37,  can  be  written  as  functions  of  t,  Xg,  Xg  and  Axs> 

,  J  '  (38) 

A  r  A  Xfj  &  /  LX} )  m 

When  there  is  no  reflection,  PQ  and  P^  are  both  large,  i.e., 

**  o(£)  ,  f>.  =  0  ({■)  . 

Significant  reflection  on  the  emitter  side  reduces  PQ  and  it  may 
indeed  attain  negative  values  for  sufficiently  strong  reflection. 


The  density  equation  (eq.  35)  with  the  boundary  conditions  PQ 
and  P^  is  a  linear  eigenvalue  problem;  its  solution  yields  A  and 
C  as  functions  of  PQ  and  The  calculated  results  are  shown  in 
fig.  2.2.2.  Since  A(t)  is  function  of  t  from  the  ionization 


kinetics,  the  value  of  t  is  thus  determined  by  a  function  of  Pg 
and  P^.  The  plasma  resistance,  R,  also  can  be  expessed  in  terms 
of  functions  of  Pg  and  B1  through  A  and  C  using  eq.27: 


~  Q  s  ioC 

f**  J  A 


(39) 


Using  the  sheath  results  which  provide  j,  Q,  PQ  and  P^,  the 
isothermal  formulation  is  complete.  The  results  are  summarized 
below.  The  quantities  Pg,  P^,  Q  and  j  are  found  from  the  sheath 
calculations  as  functions  of  t,  Xg,  Xg,  and  Axs>  i.e., 

f,.  -  f>. 

f.,- 

Q  =  Q  (T,  *-r,  LX>), 

j  C.T,  *r,%.  toQ. 
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Figure  2.2.2  The  Eigenvalue  Problem 


From  the  eigenvalue  problem  for  the  plasma  density  we  then  find 

ACr)  -  (40) 

From  the  continuity  equation  for  current  we  find 


And  from  the  electron  momentum  equation  we  find 

,  .  l(T-0  i;., 

*  Tin  (.-&KTJ  ■*■>*■'  —J . 

These  three  previous  equations  determine  X g,  *c  and  t  when  AXg  is 


(42) 
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given.  This  set  of  equations  is  valid  for  all  Axg.  Even  in  the 

case  of  Ax  <=  0  when  there  is  no  reflection,  the  calculations 
s 

differ  from  previous  isothermal  calculations  because  the 
Boltzmann  assumption  on  the  electrons  is  not  used  as  indicated  by 
the  presence  of  and  a^. 
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CHAPTER  3:  SHEATHS 


3.1  Assumptions 

3.2  General  Solution  Condition 

3.3  Equations 

3.4  Sheath  Solutions 


In  this  chapter  the  collisionless  sheath  formulation  is 
developed.  The  sheath  formulation  and  its  results  are  essential 
to  thermionic  convertor  formulation  because  two  pieces  of 
information  are  required  from  the  sheath: 

1)  the  potential  change  through  the  sheath,  and 

2)  plasma  density  boundary  conditions,  and  in  the 
non*isothermal  case,  electron  temperature 
boundary  conditions. 

Both  of  these  have  significant  effects  on  predicted  thermionic 
convertor  performance,  and  therefore  sheath  formulation  is  of 
great  interest.  In  addition  to  calculating  more  accurate  sheath 
results  than  in  previous  sheath  formulations,  we  examine  three 
exoected  emitter  sheath  phenomena: 
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1)  plasma  ions  reflected  by  the  emitter  sheath, 

2)  ions  trapped  in  the  double  emitter  sheath,  and 

3)  surface  emission  ions  from  the  emitter  surface. 

These  three  phenomena  and  more  accurate  results  require 
generalizations  to  previous  formulations  along  two  lines.  First, 
the  assumption  of  cold  ions  (  monoenergetic  plasma  ions  )  cannot 
be  retained  since  ion  reflection  begins  with  the  lowest  energy- 
plasma  ions  and  reflects  only  the  part  of  the  ion  distribution 
that  has  energies  lower  than  the  potential  rise  through  the 
sheath.  Second,  previous  formulations  have  idealized  electron 
dynamics  in  the  sheath  under  the  assumption  of  a  large  sheath 
height.  In  fact  the  actual  sheath  heights  are  of  order  unity, 
and  therefore  this  is  a  poor  assumption.  Complete  calculation  of 
the  electron  dynamics  converts  the  simple  algebraic  equations  of 
the  large  sheath  height  formulation  to  equations  containing 
integrals  over  distribution  functions  that  must  be  integrated 
numerically.  By  contrast,  the  finite  temperature  ion 
distribution  causes  theoretical  difficulties  which  must  be 
carefully  considered.  The  cold  ion  formulation  uses  the  Bohm 
criterion  on  ion  kinetic  energy  to  assure  a  self-consistant 
sheath  solution.  We  show  that  a  finite  temperature  plasma  ion 
distribution  must  have  its  low  energy  tail  cut  off  before  it 
enters  the  collisionless  sheath  since  otherwise  no  self- 
consistant  sheath  solution  exists.  We  expect  a  transitional 
region  to  exist  between  the  collisionless  sheath  and  the  neutral 
plasma  since  we  are  asymptotically  matching  the  two.  This  region 
would  be  collisior.al  but  with  relatively  strong  electric  field 
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that  would  accelerate  and  depopulate  the  ion  distribution  of  low 
velocity  members.  Therefore  we  cut  off  the  low  energy  tail  of 
the  ion  distribution.  We  show  that  the  minimum  ion  distribution 
shift  speed  depends  logarithimically  on  an  arbitrarily  chosen 
cut-off  of  low  energy  ions,  thus  we  can  therefore  pick  an  order  - 
of -magnitude  cut-off  with  little  effect  on  the  physical  outcome. 
In  both  the  cold  ion  and  finite  temperature  ion  cases,  self- 
consistantcy  is  critical  at  the  plasma-sheath  interface,  and  we 
develop  a  generalized  self-consistancy  condition  which  includes 
the  original  Bohm  criterion  as  a  special  case. 

When  trapped  ions  and  surface  emission  ions  are  considered, 
the  point  of  critical  self-consistancy  moves  away  from  the 
plasma-sheath  interface  into  the  sheath.  The  new  generalized 
self-consistancy  condition  remains  valid  for  this  case.  The 
amount  cf  surface  emission  ions  present  in  the  sheath  is 
determined  by  the  Saha -Langmuir  equation  and  thus  by  the 
temperature  and  work  function  of  the  electrode  along  with  the 
neutral  pressure  next  to  the  electrode.  On  the  other  hand,  the 
amount  of  trapped  ions  is  determined  by  the  coll-'sional  processes 
(  largely  charge  exchange  )  in  the  sheath.  The  amount  of  trapped 
ions  cannot  be  easily  calculated  because  of  the  complexity  of  the 
ccllisioi.al  processes.  Threrefore,  we  estimate  the  amount  of 
trapped  ions  and  produce  results  for  a  range  of  values. 


3 . 1  Assumptions 


-4- 


CHAPTER  3 


In  this  section  we  list  and  discuss  the  assumptions  of  the 
present  sheath  formulation  and  their  justifications.  The  present 
formulation,  which  removes  some  previous  idealizations,  assumes: 

1)  the  sheath  is  collisionless, 

2)  the  electron  distributions  from  the  emitter  and  collector 
are  Maxwellian  with  the  emitter  and  collector 
temperatures  respectively, 

3)  the  plasma  electrons  are  Maxwellian  with  the  plasma 
electron  temperature, 

4)  the  ion  distribution  from  the  plasma  has  the  neutral 
temperature,  and  is  a  shifted  Maxwellian 
distribution, 

5)  at  the  plasma-sheath  interface,  charge  neutrality  exists, 

6)  at  the  plasma-sheath  interface,  the  electric  field  is 
small  (  in  the  match-asymptotic  sense  ), 

7)  the  trapped  ions  have  a  Maxwellian  distribution  with  the 
emitter  temperature  and  an  arbitrarily  specified  density, 
and 

8)  the  surface  emission  ions  are  emitted  with  a  Maxwellian 
distribution  governed  by  the  Saha-Langmuir  equation. 

In  contrast,  most  previous  formulations  have  assumed: 

1)  the  ion  distribution  is  cold, 

2)  the  plasma  electron  density  in  the  sheath  is  Boltzmann 
with  potential,  and 

3)  the  emitter  electron  distribution  is  also  cold. 

These  three  assumptions  simplified  the  mathematics  considerably, 
but  produce  relatively  large  errors  in  the  sheath  results. 
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particularly  for  small  sheath  heights.  Since  the  emitter  sheath 
effects  of  interest  here,  namely  reflection,  trapped  ions,  and 
surface  emission  ions,  produce  smaller  sheath  heights,  it  is 
essential  to  remove  these  assumptions.  Additionally,  the  cold 
ion  assumption  is  totally  incompatible  with  ion  reflection. 
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the  continuum  plasma.  This  layer  would  also  accelerate  low 
energy  ions  and  therefore  cut  off  the  bottom  of  the  ion 
distribution.  The  cut-off  is  arbitrarily  determined,  but  has  a 
weak  logarithimic  effect  on  the  sheath  results  (  it  raises  ion 
distribution  shift  speed  ). 

Assumptions  2  and  8  regarding  the  surface  emission  of 
electrons  and  ions  are  expected  to  be  very  good  because  the 
surface  temperatures  are  high.  The  emitted  distributions  should 
be  indistinguishable  from  Maxwellian  distributions.  The  other 
assumptions  regarding  distributions,  3,4  and  7  are  more  tenuous. 
The  ion  distribution  is  assumed  to  have  the  neutral  temperature 
because  there  is  large  charge-exchange  collision  cross-section 
between  the  ions  and  the  neutrals.  The  shift  in  the  ion 
distribution  (  assumption  4  )  is  to  be  determined  as  the  minimum 
shift  required  to  construct  a  self-consistant  sheath.  The 
assumption  that  the  trapped  ions  are  at  the  neutral  temperature 
is  also  based  on  the  predominance  of  charge  exchange  collisions. 
The  assumption  that  the  plasma  electrons  are  in  a  Maxwellian 
distribution  with  their  own  temperature  is  based  on  the  length  cf 
the  convertor  plasma  region  (  approx.  15  mean  free  paths  ).  The 
convertor  is  long  enough  that  the  electrons  have  time  to 
equilibrate  with  each  other  but  not  long  enough  to  equibrate  with 
the  ions.  1 


i 


See,  for  instance,  Montgomery,  p.33. 
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3.2  The  General  Solution  Condition 

Before  proceeding  on  to  discuss  the  possible  sheath 
configurations  and  their  specific  equations,  we  develop  the 
general  solution  condition  which  applies  to  all  collisionless 
sheaths  to  preclude  consideration  of  non-self-consistant 
solutions.  Asymptotically  matching  the  collisionless  sheath  to 
the  quasi-neutral  plasma  always  requires  some  condition  on  the 
ion  distribution  function  coming  into  the  sheath  from  the  plasma. 
In  the  past  this  had  been  the  Bohm  criterion  which  assumed  cold 
plasma  ions  and  required  the  monoenergetic  plasma  ion 
distribution  to  be  shifted  up  in  velocity  to 


where  T&  is  plasma  electron  temperature,  k  is  Boltzmann's 
constant  and  K  ion  mass.  In  this  section  it  will  be  shown  that  a 
general  condition  for  solving  the  collisionless  sheath  with  a 
neutral  plasma-sheath  interface  exists  and  that  the  Bohm 
criterion  and  other  local  (  at  the  plasma-sheath  interface  ) 
matching  conditions  are  special  cases  of  the  general  condition. 
We  show  that  local  matching  is  a  necessary  but  not  sufficient 
condition  on  the  sheath  solution.  In  the  absence  of  trapped  ions 
or  surface  emission  ions,  as  is  the  case  with  most  past 
calculations,  local  matching  proves  to  be  also  sufficient. 
Further,  we  show,  in  the  case  of  no  trapped  ions  or  surface 
emission  ions,  that  for  finite  temperature  the  ion  distribution 
must  have  no  zero  velocity  ions  when  it  enters  the  collisior.less 
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A  self-consistant  sheath  solution  is  found  by  integrating  the 
nondimens ional  poisson  equation  from  the  last  section, 

dlX.  _  f/%) 

■37* ~  (1) 

The  specific  forms  of  F(x)  are  developed  in  the  next  section  but 
we  need  not  know  them  until  we  wish  to  evaluate  specific  cases. 
By  convention  X  =  0  at  the  plasma-sheath  interface  and  increasing 
X  repels  electrons.*  To  construct  a  non-trivial  solution  *  we 
integrate  eq.l  as  follows, 


J1 


Transforming  eq.3  into  eq.4  implicitly  assumes  X($)  is  monctonic 
or.  the  domain  Since  at  the  plasma-sheath  interface,  we 
assume  charge  neutrality  and  zero  electric  field,  we  have  F(C)  = 
0  and  dx/d£  =  0.  Therefore  we  can  write, 


($f=  /W 


-  +JS 


1  Since  F(0)  =  0  by  charge  neurality,  there  is  alway  the  trivial 
solution  x  =  0. 
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Equation  6  is  the  method  for  constructing  non-trivial  solutions. 
From  eq.6  and  our  previous  observation  about  converting  eq.3  into 
eq.4  we  can  see  that  that  inequality, 

CfuUx  to  ,  o  <x<Xm  m 

Jo 

must  hold  where  x  is  the  first  maximum  or  minimum  X  reaches  in 

B 

the  sheath  since  otherwise  x(£)  is  not  aonotonic.  If  we  attempt 
to  construct  a  scluiior.  which  does  not  Beet  this  condition,  then 
eq  6  causes  x£  tu  double  back  before  reaching  its  full  sheath 
height 

Snowr  be  1  x.  :r  ‘  .is.  through  3.2.1(d)  are  the 

ex:  acted  e-  s  •  ••  ••  •  *  at  icr.s  and  their  solution 

cc:.d:t  f  r  — r  «■  ,  '  arranged  in  the  order  that  we 

expect  t.'.eir  t  „*  »  •  -•  der  ;ty  through  the  convertor 

is  reduced  and  the  ra'.t  vf  erc.'ter  Richardson  current  density  to 
net  currant  density  Froir.  the  general  solution 

condition  <  eq  ,  we  can  oer.ve  r.ecesscry  local  (  at  the  plasma* 
sheath  interface  )  matching  conditions,  of  which  generalized  Bohm 
criterion  is  a  special  case.  We  can  expand  F(x)  around  x=0  as 

S.  .  . 

,  (6) 

H»t 

since  Ffx)  may  be  represented  asymptotically  in  this  form  for 
cases  we  intend  to  develop.  Then  eq.8  may  be  integrated, 


y  Qr\  y  nhr\ 

Zi+i  * 


/  /  .  t  Ihi  1  1 


(S) 
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generalized  Bohm  criterion.  The  usual  Bohm  criterion  assumes 
F(x)  is  expandable  as, 


(1C) 


in  which  case  the  criterion  is 


*’ULi0- 


(11) 


When  only  cold  ions  are  considered,  the  Bohm  criterion  (eq.ll)  is 
sufficient,  but  when  finite  temperature  ions  are  considered,  we 
must  apply  the  generalized  Bohm  criterion.  If  trapped  or  surface 
emission  ions  are  present,  we  find  that  local  matching  is  not 
sufficient  and  that  the  full  solution  condition  must  be  applied. 


3 . 3  Equations 

In  this  section  we  develop  the  sheath  equations  for  the 
collisionless  emitter  sheath  configurations  shown  in  fig.  3.2.1 
and  for  the  collisionless  collector  sheath.  In  order  to  make  the 
derivations  as  clear  as  possible,  we  divide  this  section  into 
subsections  for  each  of  the  cases  listed  previously  and  then 
follow  a  standard  format  in  presenting  the  equations.  We  divide 
each  subsection  into  the  following  order  of  presentation: 

1.  Sheath  Configuration 

2.  F . (x)  and  F  (x) 

r  e 

3.  Integrals  of  F^X)  and  Fg(x) 

4.  Equations:  F^(0)  =  1, 

r  fen  =  i . 
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niin[  [  FCQ^X  ,  0  <=  X  <=  XF  ]  =  0, 

=  o 

Jo 

5.  Useful  Resultant  Quantities 

We  begin  each  subsection  with  the  sheath  configuration  and 
appropriate  definitions,  and  then  derive  the  functions  F^Cx)  and 
Fg(X)  accordingly.  These  functions  are  respectively  the  total 
ion  and  total  electron  densities  as  functions  of  sheath 
potential.  The  the  third  part  of  each  subsection,  we  find  the 
integrals  of  F^(x)  and  Fg(x).  The  fourth  part  presents  the 
equations  which  are  solved  simultaneously  for  the  sheath.  The 
equations  F^(x)  =  1  and  Fg(x)  *=  1  represent  charge  neutrality  at 
the  plasma-sheath  interface:  the  charge  densities  there  are  both 
set  equal  to  1  because  they  are  nondimens ionalized  by  total 
plasma  density,  nQ,  at  the  plasma-sheath  interface.  The  third 
equation  is  the  general  solution  condition,  and  the  fourth 
equation  is  a  condition  which  applies  only  in  the  case  of  a 
double  sheath;  it  means  that  d$/dx  =  0  at  the  sheath  peak.  This 
can  be  seen  from  eq.  3.2.5.  These  four  equations  are  used  to 
solve  for  the  four  variables: 

N\  =  density  of  ions  moving  toward  the  emitter 
at  the  plasma-sheath  interface, 

N'k  =  density  of  electrons  crossing  the  motive 
peak  (  if  any  )  coming  from  the  emitter, 

=  density  of  electrons  moving  toward  the  emitter 
at  the  plasma-sheath  interface,  and 


u 

s 


=  shift  in  the  ion  velocity  distribution  coming 
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from  the  plasma  at  the  plasma-sheath  interface. 

The  first  three  of  these  variables,  denoted  by  capital  N's  are 
nondimensionaliztd  by  plasma  density  at  the  plasma-sheath 
interface,  ng.  In  the  cases  where  no  motive  peak  exists,  the 
fourth  equation  does  not  apply  since  N^.  is  known.  In  general 
these  equations  must  be  solved  numerically,  and  the  numerical 
methods  employed  are  straight  forward.  Therefore,  no  discussion 
of  the  methods  is  presented.  Finally,  in  the  fifth  part  of  each 
subsection,  we  present  the  useful  quantities  that  the  sheath 
results  produce. 


a.  single  electron  repelling  sheath 


1.  Configuration 


In  this  sheath,  which  occurs  at  high  current  densities,  the 

Richardson  emitted  density,  N  ,  is  equal  to  the  emitted  density, 

K 

Nj.,  because  no  motive  peak  exists, 


2.  F,(X)  and  Ffi(x) 


The  ion  density  is 
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cu-Uij*  u(/m 


u.»* 


/ « 


e(u-“^ 


(1) 


where  u  is  the  ion  distribution  shift  speed  and  u  is  the  cut- 
s  cut 

off  of  low  energy  ions.  Both  u  and  u  _  are  nondimensionalised 

s  cut 

by 


where  U  is  velocity.  The  electron  density  is 


where  .N'j.  =  and  both  are  the  nondimensional  emitted  electron 
density  at  the  plasma-sheath  interface.  The  quantity  x  is  the 
nondimens ional  plasma  electron  temperature  given  by  the  ratio 
T  /T_ .  The  first  term  in  F  (x)  represents  the  decreasing  density 
contribution  of  emitted  electrons  as  they  accelerate  down  the 
sheath.  The  second  term  represents  the  density  contribution  of 
the  plasma  electrons  which  are  decelerated  as  they  climb  the 
sheath . 


3.  Integrals  of  F.(x)  and  F  ( X ) 
i  e 
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These  are  just  the  integrals  of  eqs.  1  and  2.  The  integral  of 
F.(X)  is 


4.  Equations 


Charge  neutrality  for  the  ions,  F^(0)  =1,  is  immediately 
satified  by  eqs.  1  because  no  ions  are  reflected  and  therefore 

=1.  Charge  neutrality  for  the  electrons  is 


To  complete  the  formulation,  we  apply  local  matching  (  the  Bchm 
criterion  )  at  the  plasma-sheath  interface.  This  can  be  done 
because  1)  the  fourth  equation  does  not  apply  since  there  is  no 
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motive  peak  in  this  sheath,  and  2)  local  matching  satisfies  the 
general  solution  condition  is  this  case.  Local  matching  is 
sufficient  because  no  trapped  ions  exist  and  we  have  assumed  that 
there  is  no  surface  emission  of  ions  in  this  case.  The 
justification  for  assuming  no  surface  emission  and  carrying  out 
local  matching  in  this  case  is  that  when  this  sheath  occurs 
surface  emission  is  entirely  negligible.  1  Also,  this  case, 
because  of  simplicity,  is  used  to  explain  the  need  for  a  cutoff 
of  low  energy  ions  in  local  Bohm  matching.  Local  matching  is 
done  by  asymptotically  expanding  F(x)  =  F^x)  -  Fg(x)  at  X  =  0. 
The  lowest  order  terms  are: 


F(*)  =  o-k 


f <4  u..i  =  °7 


<6> 

»  <i 


The  X  '  term  in  this  expansion  comes  entirely  from  the  low 
energy  ions  entering  the  sheath.  It  represents  the  acceleration 

of  zero  velocity  ions  of  finite  velocity.  Since  local  matching 

requires  the  lowest  order  term  to  be  greater  than  or  equal  to 

zero  (  eq.  3.2.10  ),  the  cutoff  of  low  energy  ions,  ucut  must  be 

greater  than  zero.  Otherwise,  no  self-consistant  sheath  solution 

can  be  constructed.  We  choose  an  arbitrary  cut-off  of 

u  =  ( . l)u 
cut  mono 

(  u  is  defined  at  the  beginning  of  section  3.2  )  for  two 

mono  °  ° 

reasons:  first,  the  effect  of  the  cut-off  on  u  ,  the  ion 

s 


1  At  the  high  current  densities  corresponding  to  this  sheath, 
surface  emission  ion  density  is  10  ^  of  total  ion  density. 
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distribution  shift,  is  logarithmically  weak  (  see  fig.  3.4.3  ), 
and  second,  we  expect  an  intermediate  asymptotic  region  to  exist 
between  the  collisionless  sheath  and  the  neutral  plasma  that 
would  accelerate  the  ion  distribution  and  depopulate  it  of  low 
energy  ions . 


Vith  this  cutoff  established,  the  expansion  of  F(x)  is  then: 


where  the  term  in  brackets  is  3F./3X  and  the  generalized  Bchm 
criterion  is 


£L  - 

Equations  8  and  5  written  as, 


(6) 


/=  N'Ctoe,)+n-c,(t,tr)l 


Cs (Xtp-)  -  C(  ( X ,)  *■  BCu,),  (9) 

are  solved  numerically  for  ug  and  when  N^,  t,  and  x_  are 

given. 


5.  Useful  Resultant  Quantities 
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the  plasma-sheath  interface  moving  toward  the  emitter  to  the 
total  electron  density  there, 


(16) 


b. 

1. 


and  c.  double  sheath 


Configuration 


A 


x 


(  Ax  may  be  larger  or  smaller  than  XR  ) 


V.'e  develop  cases  b  and  c  together  because  the  differences  are 
minimal.  In  this  case  the  emitted  density,  NR,  is  no  longer  the 
same  as  the  emitted  density  crossing  the  motive  peak,  NR,  because 
the  back  sheath  height,  Ax,  repels  some  of  the  emitted  electrons. 


2.  F . (X)  and  F  (x) 
l  e 

First,  we  consider  F^x).  This  is  the  mcsr  complicated  ci  the 
two  because  we  take  into  account  three  classes  of  ions:  plasma 
ions  (  ions  coming  onto  the  sheath  from  the  plasma  ),  trapped 
ions,  and  surface  emission  ions.  The  surface  emission  ion 
density  is  known  from  the  Saha -Langmuir  equation  and  the  plasma 
ion  density  is  found  in  the  usual  way  by  setting  F^CO)  =  1. 
However,  the  trapped  ion  density  is  not  determined  by  anything  in 
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< 


I  C 


the  collisionless  sheath  theory;  it  is  determined  by  the  small 
number  of  collisions  that  do  occur  in  the  sheath.  It  is  not 
known  how  to  calculate  the  amount  of  trapped  ions  precisely,  but 
clearly  it  should  be  less  than  the  local  equilibrium  (  "fully 
trapped"  )  plasma  density.  The  major  source  of  trapped  ions  is 
expected  to  be  plasma  ions  passing  through  the  sheath  having 
charge  exchange  collisions  with  neutrals. 

To  simplify  the  equations  for  F^Cx),  we  use  a  convention  of 
comparing  trapped  and  surface  emission  ion  density  to  the  plasma 
ion  density  as  illustrated  in  fig.  3.3.1.  Figure  3.3.1  shows  the 
distribution  of  ion  velocity  at  the  sheath  motive  peak.  The 
trapped  ion  region  is  limited  by  the  lesser  of  Ax  or  Xj.  by 
definition;  an  ion  is  not  trapped  if  it  has  an  energy  greater 
than  this.  Trapped  and  surface  emission  ion  densities  are 
denoted  by  f  and  f  .  When  f  =1  and  f  =1,  the  dotted 
curve  representing  a  thermal  equilibrium  distribution  (  based  on 
the  density  of  plasma  ions  )  is  filled  by  the  trapped  and  surface 
emission  ions. 

F.(X)  is  then 


Fi(x)'  Fir  Fsirfr)t  F^tx) 


CD 


where 


FM  = 


%c 


-u} 


(2) 


W.  -V  i 


Figure  3.3.1  Trapped  and  Surface  Emission  Ions 


at  the  Sheath  Motive  Peak 
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In  these  equations  two  conventions  have  been  used  to  simplify  the 
notation: 

1)  if  a  square  root  has  a  negative  argument,  then  it  is  zero 

2)  if  an  integral's  upper  limit  is  less  than  the  lower 
limit,  then  the  integral  is  zero. 

Tnis  convention  is  followed  for  all  integrals  over  distribution 
functions.  As  before  the  densities  are  nondimens ionalized  by 
total  plasma  density  at  the  plasma-sheath  interface.  Finally  the 
parameter  fgur  is  related  to  nondimens ional  emitted  density  by 


f 

sur 


2eAX‘XE  F 


surface . 


(5) 


Electron  density,  Fe(x),  is 


The  first 

term  in  this 

equation  is 

the 

density 

of  plasma 

electrons 

in  the  sheath. 

The  second 

term 

is  the 

density  of 

emitted  electrons  in  the  sheath. 
3.  Integrals  of  F^x)  and  F^Cx) 


For  the  ions, 
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(7) 


9 


where 


l^Mdt  =  I  e™')\LL(j^_u)dli 

Utici 

/■  \jfaX~Xg 

+  |  -dM-tij)1  _ 

C^uStT-  u)  dm 

tf««t  '  ) 


(8) 
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roobi  = 

wf  + 


■for- 


(10) 


For  the  electrons, 


4.  Equations 

First,  we  have  F^O)  =  1  from  eqs.  1,2,3  and  4, 


f  e~u  diu.  4 

^  t; 


w 


(12) 


A  term  for  trapped  ions  does  not  appear  in  this  equation  because 
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trapped  ions,  by  definition,  do  not  exist  at  the  plasma-sheath 
interface.  And,  we  have  Fe(0)  =  1  from  eq.  6, 


The  other  two  equations  which  apply  in  this  section. 


rrin 


OfXi  Xg  - 


0 

) 


(14) 


and, 


FMdK=  0  (15) 


are  straight  forward  applications  of  the  integrals  of  F^(x)  and 


-27- 


CHAPTER  3 


Fe(x).  We  will  not  write  out  the  details  of  these  equations 
since  they  cannot  be  solved  analytically.  The  numerical  solution 
of  these  is  simple  in  principle:  we  have  to  solve  four 

simultaneous  equations  in  four  variables,  Nq  ,  N^,  fp2asma  and  us 

where  f  »  »  X-,  Ax,  u  _  and  t  are  given.  And,  the  four 

equation  can  be  reduced  to  two  by  using  eqs.  12  and  13  in  14  and 

15  to  eliminate  f  ,  and  N_  .  However,  the  actual  solution 

plasma  0 

for  each  case  requires  approximately  15  sec  of  CPU  time  on  an  IBM 
370-3081. 

5.  Useful  Resultant  Quantities 

The  results  of  this  sheath  are  presented  in  the  same  way  as  in 

3.3(a)5  except  that  Nj.  is  a  dependent  variable  and  the  results 

now  depend  on  f .  ,  f  and  Ax : 

K  tr  sur 

U  =  u  (1*J  AX,  T,  fir,  -hr); 

j  *  j  ( ^e)tX)T,  ftn4atr) j 

Q-  Q  (  jT) ‘flrj-tsur), 

=  Vote, MtTjir,  hr), 


(16) 
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In  this  case,  we  have  reversed  the  sense  of  X  for  convenience. 


2.  Fi(X)  and  Fg(x) 

The  ion  density,  F^(x),  is 


We  have  contracted  this  expression  for  ion  density  to  eliminate 
hT  by  using  F^O)  =  l  immediately .  For  this  sheath 
configuration,  there  cannot  be  trapped  ions. 

The  electron  density,  Fe(X),  is 


where 


F  (6)  =  N~rhf£  =  / 


(3) 
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3.  Integrals  of  F.(x)  and  F  (X) 
1  " 


'[jrf  -  u)  du  J  m 

*  *•[£(«*  fr*'*-#)  i-P-o] 


A.  Equations 

Instead  of  presenting  only  the  equations,  in  this  case  we  can 
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produce  an  analytical  result.  We  find,  analytically,  that  the 

minimum  for  u  is  always  less  than  zero  for  small  Ax  under 
s 

certain  easily  satisfied  conditions.  Regardless  of  Ax,  numerical 
calculation  always  supports  ug  =  0. 

To  demonstrate  that  ug  =  0  for  small  Ax  we  assume  fgur  =  0  for 

simplicity  *  and  take  the  asymptotic  expansion  of  the  integrals 

for  u  =0, 
s 


The  general  solution  condition  then  requires 


No 


< 


Vz. 

t+  r'/z 


(8) 


*  In  the  cases  that  this  sheath  occurs,  f  _  is  always  small  ( 

sur 

10*  )  and  does  not  affect  the  numerical  conclusion  that  u  =  0. 


-31- 


CHAPTER  3 


We  now  show  that  for  J  >=  0  (  positive  net  current  density  )  eq. 
8  is  satified: 


(l-  K(nff'))' 


(9) 


Then  J  >=  0  implies 

N°~  -  •  (10) 

Therefore  eq.  8  is  satisfied  for  J  >=  0  and  u  =  0  is  a  self- 

s 

consistant  solution. 


5.  Useful  Resultant  Quantities 


The  resultant  quantities  of  this  sheath  must  be  presented  to 
the  rest  the  thermionic  convertor  formulation  differently  because 
there  is  no  sheath  height,  X^..  In  fact  from  eqs .  3  and  9  ( 
charge  neutrality  and  net  current  )  we  can  immediately  derive 


(id 


The  definitions  of  j  and  Q  are  as  in  3.3(a)  resultant  quantities. 
Because  of  this,  we  present  these  sheath  results  in  this  case  as 


Q-Q(j,x). 
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( 


e 


■  < 


a 


e.  collector  sheath 

1.  Configuration 


Since  the  collector  is  assumed  to  emit  nothing,  this  case  is 
particularly  simple.  As  in  all  other  cases,  the  electron 
distribution  function  is  assumed  to  have  no  shift.  However,  as 
will  be  seen  in  section  4.2  the  collector  sheath  height  is 
usually  small  and  goes  to  zero  with  ion  reflection  at  the 
emitter.  Because  of  this  it  may  be  of  interest  in  further  work 
to  allow  the  electron  distribution  to  have  a  shift  and  to 
simultaneously  allow  an  ion  repelling  collector  sheath. 


2.  F. (X)  and  F  (X) 
l  e 


The  ion  density,  F^Cx),  is 


(1) 


» 


t  J 


and  the  electron  density,  F  (X),  is 


5.  Useful  Resultant  Quantities 
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The  results  of  this  sheath  are  presented  as: 

u  =  u  (XCjr), 
or,  =  o',  (tipO , 


where  is  the  fraction  of  electrons  at  the  plasma-collector 
sheath  interface  moving  toward  the  collector, 


(7) 


3.4  Sheath  Solutions 

Again,  we  divide  this  section  into  subsections  for  each  type 
of  sheath.  Most  of  the  results  in  this  section  are  on  the  double 
emitter  sheath  because  its  is  the  most  complex  case  and  contains 
the  emitter  sheath  phenomena  of  interest  -  trapped  ions  and 
surface  emission  ions. 

a.  single  electron  repelling  sheath 

This  sheath  occurs  at  the  emitter  when  the  net  current 
density,  J,  is  greater  than  appoximately  75%  of  J^.  We  do  not 
present  results  for  this  case  because  they  are  expected  to  be 
similar  to  conventional  theories. 


b.  and  c.  double  sheath 
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This  subsection  on  the  double  emitter  sheath  contains  most  of 
the  results  of  interest.  In  this  subsection,  we  first 
distinguish  the  results  of  the  present  full  collisionless  sheath 
theory  from  past  approximate  theories  but  still  assuming  no 
trapped  or  emitted  ions .  We  then  demonstrate  the  effect  of 
trapped  ions  on  the  sheath.  We  make  no  mention  of  surface 
emission  ions  in  this  subsection  because  of  their  small  amount 
and  small  effect  in  the  convertor  at  the  current  densities  we 
make  calculations  for.  However,  surface  emission  is  fully 
included  in  all  sheath  results  used  in  the  thermionic  convertor 
calculations . 

The  essential  diferences  in  results  of  the  present  full 
collisionless  theory  and  past  approximate  theories  is  illustrated 
in  figs.  3.4.1,  3.4.2  and  3.4.3.  Figure  3.4.1  demonstrates  the 
relationship  between  sheath  height,  x^.,  and  normalized  current, 
j.  The  first  curve  (  to  the  right  )  is  the  simplest  approximate, 
assuming  that  1)  the  emitted  electrons  are  monoenergetic  at  their 
average  speed,  2)  the  plasma  electrons  are  in  a  fully  Boltzmann 
distribution  despite  a  finite  sheath  height,  and  3)  the  plasma 
ions  are  monoenergetic.  The  second  curve  removes  the  Boltzmann 
plasma  electron  assumption.  The  third  curve  removes  the  cold 
emitter  electron  assumption.  And,  the  fourth  curve  removes  the 
final  assumption  of  cold  plasma  ions.  The  first  two  changes  are 
unequivocal  corrections  to  the  sheath  theory.  However,  the  last 
correction,  removal  of  the  cold  ion  approximation  requires  the 
imposition  of  the  cut-off  of  low  energy  ions.  The  cut-off,  as 
discussed  previously,  is  set  at  10%  of  the  monoenergetic  Bohm 
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Figure  3.4.2  Shift  Speed  and  Sheath  Approximations 


Trapped  ions  have  a  substantial  effect  on  the  sheath  solution. 

Figure  3.4.4  shows  the  general  solution  condition  for  a  double 

sheath  with  xr  and  for  f  =0.0  and  for  f„  =0.2.  In  the  no 
E  tr  tr 

trapped  ion  case  (  f  =  0.0  )  we  have  local  matching  in  which 
the  critical  point  occours  at  the  plasma  sheath  interface  (  X  =  0 
)  As  trapped  ions  are  added,  the  critical  point  moves  into  the 
sheath  (  f  =0.2  ).  Figure  3.4.5  shows  the  sheath  potential 
versus  nondimens iona]  position,  £  =  x/X^.  The  two  cases  shown 
are  for  f  =  0.0  and  -  0.2  corresponding  to  fig.  3.4.4.  In 
the  no  trapped  ion  case  the  potential  drops  rapidly  from  its  peak 
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Figure  3.4.3  Effect  of  Cut-Off  on  Shift  Speed 


at  Xj.  =  2.0  to  Xg  =  0.0.  In  the  f  -  0.2  trapped  case  the 
sheath  drops  rapidly  to  the  critical  catching  point,  where  it 
levels  out  asymptotical ly  and  then  continues  its  drop  to  X  =  C.O. 
This  asymptotic  region  has  a  finite  slope  in  fig.  3.4.3  because  a 
numerical  factor  has  been  added  for  display  purposes.  The 
existence  of  this  asymptotic  ''flat  c-pot"  is  not  troublesome 
because  a  slight  increase  in  u&  will  remove  it.  Figure  3.4.6 
shows  the  trapped  ion  effect  on  u  for  various  Xr.  Only  a 
certain  amount  of  trapped  ions  can  be  added  before  no  sheath 
solution  exists.  Trapped  ions  added  for  a  given  sheath  height, 


Figure  3.4.4  Solution  Condition  with  Trapped  Ions 


Xj. ,  require  a  larger  (  emitted  electrons  )  and  a  smaller  ( 

plasma  electrons  ).  The  solution  fails  when  the  required  is 

less  than  zero  (  a  negative  density  of  plasma  electrons  ).  The 
failure  point,  in  f  ,  decreases  as  Xj.  increases.  This  is 
because  the  total  density  of  trapped  ions  that  corresponds  to  a 
given  f  rises  exponentially  with  Xv.  Figure  3.4.7  shows  the 
effect  of  trapped  ions  on  nondimens ional  current  density,  j.  The 
essentia]  feature  here  is  that  trapped  ions  reduce  X^.  for  a  given 
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Figure  3.4.6  Shift  Speed  with  Trapped  Ions 


zero  the  matching  remains  local  while  trapped  ions  accumulate  in 
sheath  peak.  This  causes  ug  to  rise  and  also  causes  j  to  rise  ( 
fig.  3.4.9  ).  For  display  purposes  ug  has  been  limited  to  3  by 
the  numerical  routine.  As  can  be  seen  ug  increases  exponentially 
with  Ax  until  Ax  =  X^  where  trapped  ions  begin  to  control 
matching.  For  the  purpose  of  clarity,  u  is  shown  to  have  a 
finite  slope  at  Ax  =  Xj.  when  in  fact  the  slope  is  infinite. 

After  AX  =  Xr  is  passed  u  drops  to  its  value  shown  in  fig. 
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Figure  3.4.8  Effect  of  Back  Sheath  Height  on  Shift  Speed 


e .  collector  sheath 

This  is  the  classical  Bohm  matching  case.  Figure  3.4.10  shows 

the  collector  shift  speed,  us>  versus  collector  sheath  height, 

Xj,.  The  curves  shown  are  for  various  ion  distribution  cut-offs, 

u  =  . 1,  .01,  .001,  .0001.  The  notable  feature  here  is  that  u 
cut  s 

drops  to  negative  infinity  at  a  finite  but  small  xr.  This  is 

w 

expected  because  a  small  sheath  should  not  have  a  pre-sheath 
region  capable  of  accelerating  ions  to  a  high  speed. 
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CHAPTER  4: 


ISOTHERMAL  SOLUTIONS 


4.1  Effects  of  Ion  Reflection 

4.2  Effects  of  Trapped  Ions 

4.3  Effects  of  Surface  Emision 

4.4  Comparison  with  Experimental  Work 
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In  this  chapter  we  develop  and  discuss  isothermal  solutions 
for  the  thermionic  convertor  with  the  emitter  sheath  phenomena  of 
ion  reflection,  trapped  ions  and  surface  emission  ions  included. 
All  three  of  these  phenomena  increase  in  significance  as  net 
current  density  through  the  convertor  is  reduced.  Each  of  the 
these  reduces  the  net  ion  loss  rate  to  the  emitter  and 
consequently  increases  arc-drop  (  therefore  degrading  performance 
at  low  current  densities  ).  This  increase  in  arc-drop  is  in 
agreement  with  the  same  tendency  in  the  experimental  results. 
However,  the  experimental  results  also  show  a  plateau  (  of  low 
arc-drop  )  at  low  current  density.  This  plateau  occurs  at  a 
current  density  corresponding  to  significant  surface  ion  emission 
and  is  therefore  thought  to  occur  as  surface  emission  replaces 
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volume  ionization  as  the  dominant  source  of  plasma  ions. 
Unfortunately,  the  theoretical  calculations  cannot  be  carried 
into  this  region  because  the  collisionless  collector  sheath 
matching  (  to  the  neutral  plasma  )  fails. 

To  provide  a  realistic  framework  for  presenting  the  results  of 
this  chapter,  we  consider  the  convertor  conditions  shown  as  case 


1 « 

CASE  1  CASE  2 

•  4 

• 

T£  =  1500  K  T£  =  1750  K 

Tc  =  750  K  Tc  =  750  K 

; 

K 

p  =1  torr  p  =1  torr 

cs  *cs 

.  j 

• 

d  =  10  mil  d  =  10  mil 

=  2.12  eV  =  2.67  eV 

-  .*• 

!• 

*c  =  1.60  eV  dc  *  1.73  eV 
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2  2 
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Figure  4.1  Isothermal  Solution  Conditions 

1  in  fig.  4.1.  Case  2  is  shown  because  it  has  the  largest 

i 

surface  emission  of  any  typical  thermionic  convertor  operating 

\S‘. 

condition  (  because  the  work  function  is  high  and  the  temperature 

j< 

is  also  high  ).  Instead  of  presenting  case  2  seperately,  we 

1 

demonstrate  the  effects  of  surface  emission  in  case  1  by 

increasing  the  surface  emission  by  a  factor  of  100  thereby 

•L 

k 

bringing  it  up  to  the  level  in  case  2.  The  net  current  density 

1 

at  which  surface  emission  becomes  significant  can  be  estimated  by 

multiplying  J^s+  by  the  square  root  of  the  ion  to  electron  mass 

i< 

1  .. 
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ratio  (  approximately  500  ) .  In  case  1  this  means  that  surface 

2 

emission  becomes  significant  at  J  =  .01  amps/cm  while  in  case  2 

2 

significant  surface  emission  begins  at  J  =  1.0  amps/cm  . 

4.1  Effects  of  Ion  Reflection 

In  this  section  we  develop  the  isothermal  results  for  case  1 
with  ion  reflection,  but  without  trapped  ions  and  with  the  small 
amount  of  surface  emission  ions  of  case  1.  Figure  4.1.1  is  the 
C-V  diagram  for  this  case.  The  dotted  line  extending  upward  from 
point  A  is  the  single  electron  repelling  emitter  sheath  solution. 
However,  we  have  nor  taken  recombination  into  account  in  this 
isothermal  calculation  nor  have  we  included  the  Schottky  effect, 
both  of  which  are  expected  to  become  important  at  current 
densities  near  JR.  Therefore  the  curve  above  point  A  should  be 
treated  as  qualitative  and  not  quantitative. 

The  interest  of  this  thesis  begins  at  point  A,  where  the 
single  sheath  doubles  over.  Between  points  A  and  B,  where  the 
back  sheath  height,  Ax,  is  less  than  the  sheath  height,  Xj.,  the 
emitter  sheath  is  non-reflecting.  In  this  region  the  sheath 
heights,  Xj.  and  X^,  remain  constant  while  the  plasma  density  is 
proportional  to  net  current,  J  (  the  normalized  plasma  density 
nc/J  is  constant  ).  Only  the  back  sheath  height,  Ax,  changes  and 
the  C-V  curve  in  this  region  is  Boltzmann  (  the  arc-drop  is 
constant  ).  Beginning  at  point  B  and  continuing  to  point  C,  the 
double  emitter  sheath  reflects  plasma  ions  because  the  back 
sheath  is  larger  than  the  front  sheath,  in  other  words  the 


CvelB) 

Figure  4.1.1  C-V  Diagram  with  Ion  Reflection  for  Case  1 


reflective  potential,  Ax  =  Ax  -  x_,  is  positive.  The  result  is 
that  net  ion  loss  rate,  u,  decreases  and  that  arc-drop  increases. 
The  dotted  curve  BD  is  the  same  double  sheath  except  that  it 
assumes  no  ions  are  reflected;  therefore  TT  is  constant  and  arc- 
drop  is  constant.  The  two  curves  BC  and  BD  are  almost 
indistinguishable  because  the  increase  in  arc-drop  is  small  until 
the  net  current  density  is  extremely  small.  The  reason  for  this 
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is  that  the  shift  speed  is  approximately  u  -2  and  therefore  a 

s 

large  increase  in  reflective  potential  is  required  to  change  u 
significantly  (  the  half  reflection  point  is  Axg  =  4.0  or 

-4  2 

approximately  J  =  J^e  =  .4  amp/cm  ). 

The  curve  EF  is  the  single  electron  repelling  emitter  sheath 

case.  It  is  the  limiting  case  for  large  amounts  of  trapped  ions 

in  which  the  double  sheath  peak  has  been  completely  suppressed 

by  the  trapped  ions.  As  is  explained  in  sections  3.3(d)  and 

3.4(d),  the  solution  condition  is  satisfied  by  u  =0.  This 

s 

curve  is  not  topologically  connected  to  the  curve  ABC;  it  will  be 
shown  in  section  4.2  that  trapped  ions  move  ABC  toward  the  single 
ion  repelling  sheath  case.  Curve  is  much  steeper  (  a  faster 
increase  in  arc-drop  )  in  this  case  because  ug  *  0  (  the  half 

2 

point  in  ion  reflection  is  approximately  J  =  8  amp/cm  ).  Curve 

EG  is  the  single  ion  repelling  case  assuming  no  reflection  and  is 

therefore  a  Boltzmann  line  with  constant  arc-drop. 

At  points  F  and  C  the  solutions  fail  at  the  collector.  The 
explanation  for  this  failure  is  best  given  by  examining  figs. 
4.1.2,  4.1.3,  4.1.4  and  4.1.5.  Figure  4.1.2  is  the  normalized 
plasma  density  through  the  convertor  gap.  The  highest  curve  with 
no  reflection,  AXg  =  0,  has  the  largest  plasma  density  at  the 
collector  but  the  lowest  plasma  density  at  the  emitter.  Ion 
reflection,  which  decreases  the  ion  loss  rate  to  the  emitter, 
raises  the  plasma  density  at  the  emitter  but  lowers  the  plasma 
density  at  the  collector.  The  lower  plasma  density  at  the 
collector  forces  a  smaller  collector  sheath  height  to  pass  the 


Figure  4.1.7  Normalized  Plasma  Density  with  Reflection 


net  current  density.  This  can  be  seen  from  eq.  2.2.10.  Figure 
4.1.3  is  the  potential  through  the  convertor  under  the  same 
reflection  conditions  as  in  fig.  4.1.2.  In  fig.  4.1.3  the  first 
two  spaces  on  the  left  make  up  the  double  emitter  sheath,  and  the 
last  space  on  the  right  is  the  collector  sheath.  The  region 
between  the  two  sheaths  is  the  neutral  plasma  region.  In  the  no 
reflection  case,  it  can  be  seen  that  the  potential  has  a 
pronounced  well  in  the  middle.  This  is  the  result  of  the  large 
plasma  density  in  the  middle.  As  reflection  increases,  this  well 
disappears  on  the  collector  side  of  the  plasma  because  resistive 


Figure  4. 1.4  Collector  Sheath  Failure 


failure  occurs.  Collector  sheath  height,  xr  goes  toward  zero, 

V 

the  shift  speed,  ugc,  goes  toward  negative  infinity  (  see  fig. 

3.4.10  ),  and  the  ion  loss  rate  to  the  collector,  IT  ,  is  driven 

c 

to  zero.  Figure  4.1.5  shows  the  changes  in  the  emitter  sheath 
height,  ion  shift  speed  and  ion  loss  rate.  When  the  collector 
sheath  failure  occurs,  the  ion  loss  rate  to  the  collector  is  zero 


(  u  *  0  )  and  the  corresponding  plasma  ion  distribution  at  the 
collector  is  bunched  at  zero  velocity  (  ugc  *  -»  ) .  While  the 
mathematics  hold  self-consistantly  until  "uc  =  0,  the  physics  is 
clearly  poor  at  this  point  because  TT.  =  0  demands  that  the  plasma 
ions  at  the  collector  have  zero  energy  (  zero  temperature  and 
zero  mean  velocity  ) .  An  estimate  of  when  the  physics  becomes 
poor  is  ugc  =  0.  At  this  point  the  net  ion  loss  rate  is  close 
to  the  thermal  speed.  A  second  physical  difficulty  that  occurs 
with  collector  sheath  failure  is  that  the  electron  Mach  number 
there,  Q  ,  (  from  eq.  2.2.10  )  becomes 
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because  the  collector  sheath  height  approaches  zero  (  actually 
about  .001  ).  In  the  present  continuum  formulation  of  the  plasma 
region,  it  was  assumed  in  eq.  2.1.13  that  Q£  is  small  so  that  the 
electron  momentum  term,  ugdue/dx  can  be  neglected. 

One  could  take  the  solution  below  the  collector  sheath  failure 

point  if  "u"  could  attain  negative  values  or  if  Q  could  attain 
c  c 

values  larger  than  A2/ir).  There  is  no  physical  basis  for 
assuming  that  uc  can  become  negative  since  the  collector  emits 
nothing.  However,  there  is  a  physical  basis  for  allotting  to 
be  larger  than  /(2/ir),  (  an  electron  distribution  shift  )  as  can 
be  seen  in  fig.  4.1.3:  the  potential  drop  nearing  the  collector 
becomes  progressively  more  electron  accelerating  as  the  collector 
sheath  fails  and  therefore  the  electron  distribution  should  De 
shifted  as  the  ion  distribution  is  in  an  electron  repelling 
sheath.  However,  this  would  clearly  invalidate  the  assumption 
that  the  electron  momentum  term  is  negligible.  Therefore  the 
momentum  term  must  be  added  to  explore  further  in  this  direction 
and  this  has  not  been  done  because  of  the  resulting  complexity  in 
the  equations. 

Comparison  of  fig.  4.1.4.  to  fig.  4.1.5  at  the  collector 
sheath  failure  point  (  Axg  =  2.5,  iT  =  0  )  shows  that  the  ion 
loss  rate  to  the  emitter  is  positive.  At  this  point  the  plasma 
is  still  ignited  and  generating  ions  as  can  be  seen  from  figs. 
4.1.6  and  4.1.7.  The  ionization  coefficient,  A,  has  dropped  by 
50%,  but  the  plasma  electron  temperature  has  dropped  by  only  5%. 
Finally,  we  note  in  fig.  4.1.6  that  the  normalized  plasma 


» 
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This  can  be  seen  from  eq.l.  However,  as  we  see,  the  collector 
failure  occurs  before  t  has  dropped  more  than  5%.  Consequently 
we  do  not  see  any  plateau  or  decrease  in  arc-drop  as  net  current 
density  is  decreased  in  the  present  calculations. 


definition  in  section  3.3.  )  Curve  AHIJ  is  the  C-V  characteristic 


for  f  =  0.10  At  point  A  there  cannot  be  any  trapped  ions  since 
the  back  sheath  height.  Ax,  is  zero.  Therefore  the  trapped  C-V 
merges  into  the  non-trapped  curve  there.  The  actual  amount  of 
trapped  ions  on  the  f  =0.10  curve  increases  from  zero  at  point 
A  to  the  full  10%  of  a  thermal  distribution  at  point  H  where  the 
back  sheath  height,  Ax,  is  equal  to  the  sheath  height,  Xj..  The 
shift  speed  increases  on  AH  from  1.95  to  3.00.  This  correspondes 
to  what  is  seen  in  fig.  3.4.8  where  Ax  <  Xj..  The  rise  in  shift 
speed  has  been  limited  to  3.00  as  in  fig.  3.4.8.  This  limit  is 
placed  on  the  shift  speed  because  a  sheath  with  height  of  about 
1.0  should  not  have  a  pre-sheath  region  capable  of  shifting  the 
entire  distribution  so  far.  In  fact  limiting  the  shift  speed  is 
equivalent  to  increasing  the  cut-off  speed  for  the  ion 
distribution  function.  The  arc-drop  decreases  as  result  of  the 
increase  in  ug  and  the  consequent  increase  in  the  net  ion  loss 
rate  to  the  emitter.  A  "hump"  can  be  seen  on  AH  where  the  shift 
speed  hits  3.00.  The  arc-drop  is  lowest  on  this  "hump"  because 
the  shift  speed  is  at  its  maximum  of  3.00.  Between  points  H  and 
I  the  back  sheath  height  remains  equal  to  the  sheath  height,  Ax  - 
Xr  =  X  =0.  On  this  segment,  u  decreases  to  1.25,  therefore 
increasing  arc-drop.  The  drop  in  ug  at  Axg  =  0  can  be  seen  in 
fig.  3.4.8  also.  From  point  I  to  point  J,  the  shift  speed 
remains  constant  at  1.25  and  the  ion  loss  rate  decreases  because 
of  reflection.  The  other  trapped  cases,  f  =  0.2,  0.3,  and  0.4 
have  not  been  connected  because  they  hit  the  3.00  maximum  shift 
speed  much  sooner  than  in  the  f  =0.1  case  as  can  also  be  seen 
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from  fig.  3.4.8.  Point  J  is  the  collector  sheath  ^ai^ure  point. 
Each  of  the  f^  =0.2,  0.3  and  0.4  curves  begins  -  ^ 

ends  at  the  collector  sheath  failure  point.  It  should  he  noted 
that  each  of  the  trapped  ion  curves  failes  at  a  higher  current 
than  the  last  because  the  shift  speed  is  lower. 


4.3  Effects  of  Surface  Emission 


In  this  section  the  effect  of  surface  emission  is  discusset*- 
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Figure  4.3.1  C-V  Diagram  with  Surface  Emission 
and  Trapped  Ions 
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Figure  4.3.1  adds  the  effect  of  surface  emission  to  fig.  4.2.1. 
On  the  f  =  0.10  curve,  surface  emission  is  added  by  multiplying 
the  actual  small  amount  of  surface  emission  in  case  1  by  a  factor 
of  100.  This  brings  the  surface  emission  up  to  the  level  in  case 

2 

2,  making  it  significant  at  J  =  1.0  amp/cm  .  It  can  be  seen  that 
surface  emission  increases  arc-drop;  it  does  so  in  exactly  the 

same  way  as  reflection  or  trapped  ions  do  -  it  decreases  the  net 

loss  rate  of  ions  to  the  emitter.  Also  the  collector  sheath 

failure  occurs  at  point  K  in  exactly  the  same  way  as  in  section 

4.1. 

4.4  Comparison  with  Experimental  Work 

Figure  4.4.1  puts  the  isothermal  results  of  fig.  4.3.1  next  to 
the  experimental  results  of  fig.  1.1.2.  The  point  of  this 
comparison  is  that  the  steepness  of  the  C-V  characteristic  in  the 
experimental  convertors  can  be  explained  by  a  decreasing  ion  loss 
rate  to  the  emitter.  We  have  shown  that  all  three  of  the 
expected  emitter  sheath  phenomena  decrease  the  ion  loss  rate  to 
the  emitter.  We  cannot  calculate  the  amount  of  trapped  ions  in  a 
collisionless  sheath  without  knowledge  of  the  collisional 
processes.  However,  the  experimental  C-V  suggests  that  if  the 
amount  of  trapped  ions  (  f  )  increases  from  0%  at  J  =  14 
2 

amps/cm  (  the  double  sheath  formation  point  )  to  30%  at  J  “  2 
2 

amps/cm  (  the  collector  sheath  failure  at  f  =0.30  )  then  the 

steepness  could  result  from  trapped  ions  reducing  the  ion  loss 
rate  to  the  emitter.  Since  these  percentages  are  based  on  a 
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Figure  4.4.1  Isothermal  versus  Experimental  C-V  Diagrams 


emitter.  Since  these  percentages  are  based  on  a  thermal 
distribution  of  ions,  they  seem  physically  reasonable. 
Unfortunately,  the  collector  sheath  failure  prevents  us  from 
going  to  the  point  ir.  the  calculations  where  t  drops  enough  to 
make  surface  emission  the  source  of  ions. 

The  experimental  curve  is  nearly  a  constant  .05  volts  below 
the  isothermal  result  (  f  =0.10  )  except  at  high  current 
densities  and  at  the  "hump".  Comparison  of  the  curves  at  high 
current  density  is  not  valid  since  neither  the  Schottky  effect 


effect  nor  recombination  have  been  included.  The  Schottky  effect 

2 

is  important  above  12  amps/cm  in  this  case  because  the  emitter 
sheath  is  single  electron  repelling  (  to  the  plasma  )  and 

therefore  puts  a  strong  electric  field  against  the  emitter  with 

the  appropriate  sign.  Recombination  is  also  potentially 

important  because  the  plasma  density  scales  with  current  density 

and  at  high  current  densities  the  plasma  density  in  the  middle  of 

the  convertor  approaches  the  Saha  density.  The  .05  volt 

difference  may  or  may  not  be  explained  by  a  discrepency  in  the 

assumed  collector  work  function.  At  800  K  the  collector  emits 

essentially  nothing  and  therefore  any  change  in  the  collector 

work  function  directly  affects  output  voltage.  If  the  collector 

work  function  were  in  fact  1.65  vclts  instead  of  1.60  volts  then 

the  isothermal  result  would  lie  nearly  on  top  of  the  experimental 

result.  We  have  not  adjusted  the  assumed  collector  work  function 

so  as  to  illustrate  the  importance  it  and  therefore  the 

importance  cf  the  surface  physics  of  the  adsorbed  cesium  layer. 

The  "hump"  should  not  be  taken  as  an  expected  experimental  result 

since  it  results  from  the  interaction  of  the  trapped  ions  with 

the  plasma-emitter  sheath  interface  (  fig.  3.4.8  ).  Instead  it 

should  be  taken  as  a  second  reason  (  in  addition  to  the  cut-off 

of  the  ion  distribution  )  for  further  study  of  the  matching 

region  between  the  collisionless  sheath  and  the  neutral  plasma. 
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CHAPTER  5:  NON -ISOTHERMAL  SOLUTIONS 


5.1  The  Implicit  Computational  Scheme 

5.2  Non-Isothermal  Results 


The  non-isothermal  solution  for  the  thermionic  convertor  is 
found  by  integrating  (  through  time  )  the  set  of  parabolic 
Partial  Differential  equations  (  eqs.  2.1.14  and  2.1.15  ): 


These  equations  can  be  marched  forward  in  time  with  an  explicit 
scheme  by  computing  all  of  the  quantities  on  the  right  at 
previous  time  step  (Lawless).  This,  however,  encounters  a 
stability  limit  on  At  which  slows  the  numerical  solution 
excessively.  The  energy  equation  (  eq.  2  )  has  a  time  scale 
inherently  500  times  faster  than  the  diffusion  equation  (  eq.  1  ) 
because  the  mass  ratio  of  electron  to  ions  is  approximately  500. 
The  Stability  limit  for  the  faster  energy  equation: 
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(3) 


where  Ax  is  the  space  step  size,  te^ec  is  the  characteristic 

electron  energy  time  scale  and  Atgtaj}  is  the  stability  limit  time 

step  size.  Since  the  energy  equation  is  the  faster  of  the  two  in 

the  system,  its  stabilty  limit  controls  the  time  step  size  for 

the  set  of  both  equations.  It  can  be  seen  that  if  10  space  steps 

are  used,  then  At^^  =  100t  j  .  Since  we  want  to  find  a  steady 

state  solution,  we  need  about  10t..,,,  but  t, =  500t  ,  ( 

diff  diff  elec 

approximately  ),  therefore  we  would  need  about  5  x  10^  time  steps 
to  achieve  a  steady  state  solution. 


The  implicit  scheme  does  not  have  the  stability  limit  and 

3 

therefore  needs  only  5  x  10  time  steps  to  achieve  the  steady 

4 

state.  Each  time  step  requires  about  10  floating  point 
operations  (  flops  )  whether  or  not  the  implicit  scheme  is  used  ( 

assuming  10  space  steps  ).  For  the  explicit  scheme  this  would 

9 

result  in  5  x  10  flops  while  the  implicit  scheme  would  require 
7  4 

only  5  x  10  flops.  In  both  cases  the  estimate  of  10  flops  in 
each  time  step  assumes  that  the  sheath  calculations  are  a 

negligible  part  of  this.  In  fact  the  general  sheath  solution 

g 

requires  approximately  10  flops  in  its  full  form.  This  would 

result  in  5  x  10^  flops.  Therefore,  running  even  the  implicit 
non -isothermal  scheme  requires  that  a  sheath  model  be  used. 

Because  of  this,  exploring  thermionic  convertor  sheath  effects 

with  the  non-isothermal  theory  is  difficult  and  we  demonstrate 

only  the  simple  non-reflecting,  non-trapped  ion,  non-surface 
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emission  case  in  this  chapter  and  compare  it  to  the  isothermal 
result. 


5.1  The  Implicit  Computational  Scheme 


The  diffusion  equation  can  be  written  as 


ifi_  .218  .  L?n  .  r  j 
m.  -  a  *  oxr  +  C  -fa 


“  9*' 

and  the  energy  equation  can  be  written  as 


(1) 


S'* 


(2) 


The  quantities  a  through  h  are  not  constant  but  approximately  so. 
The  explicit  scheme  is  constructed  as  follows: 


a-iciski  +  b.dLjL+  tni+d 


Ltz 


(3) 


and 


'rf-V/ 


#  • 

=  £  ^  ’Tht~l£l  +  Q'tf  +  h 

^  A**-  2AX  j 


(A) 


At  w  24X 

where  j  is  the  time  index  and  i  is  the  space  index.  If  a  and  e 
are  the  only  non-zero  coefficients  (  a  and  e  dominate  ),  then  the 
stability  criteria  are 


.  A*1 


(5) 


and 


&4  <  J  A*1 
z  e 


(6) 


We  will  not  analyse  the  stability  criteria  in  any  more  detail, 
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but  simply  develop  the  implicit  scheme  which  has  no  stability 
limits.  The  implicit  scheme  is  constructed  by  taking  the 
variables  at  the  advanced  time  step  in  the  right  sides  of  eqs.  1 


and  2: 


=  q£'?j r^tb^  *  cnr-td 


These  equations  can  be  written  as 


t*&*  £  tt-  fe-ggfc1 

(10) 

which  form  tridiagonal  matrices,  that  with  the  boundary 
conditions,  can  be  inverted  in  3N  operations  where  N  is  the 
number  of  space  steps . 


5.2  Non-Isothermal  Results 


Figure  5.2.1  shows  the  electron  temperature  distribution 
across  the  convertor  gap  for  two  cases:  the  non-isothermal  case 
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Figure  5.2.1  Non- Isothermal  Electron  Temperature 


and  the  corresponding  isothermal  case  at  Ax  =  0  and  J  *  12 

s 

2 

amps/cm  .  The  electron  temperature  remains  close  to  2  in  the 
non- isothermal  case.  Figure  5.2.2  is  the  corresponding 

normalized  plasma  density.  It  can  be  seen  that  the  plasma 

density  in  the  isothermal  case  corresponds  closely  with -the  non- 

isothermal  case.  The  seemingly  large  difference  in  plasma 

density  in  the  middle  of  the  gap  not  important  since  its 

contribution  to  total  plasma  resistance  is  small  and  all  other 

effects  of  plasma  density  result  from  the  plasma  density  at  the 


emitter  and  collector. 


M-l50TtfEWhL 
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CHAPTER  6:  CONCLUSIONS  AND  RECOMMENDATIONS 

6.1  Conclusions 

6.2  Recommendations  for  Further  Work 


6 . 1  Conclusions 

Using  only  a  simple  isothermal  model,  the  emitter  sheath 
effects  of  trapped  ions,  ion  reflection  and  surface  emission  ions 
can  explain  the  steep  C-V  characteristic  of  the  thermionic 
convertor.  Contrary  to  intuition,  all  three  of  these  effects 
increase  arc-drop  in  the  thermionic  convertor  because  they 
increase  plasma  density  at  the  emitter  and  increase  electron 
energy  loss  to  the  emitter  to  a  greater  degree  than  they  decrease 
ionization  energy  loss.  The  low  current  density  plateau  observed 
in  thermionic  convertors  may  well  still  be  explained  by  surface 
emission,  however  we  cannot  demonstate  this  because  of  our 
collector  sheath  failure. 


6.2  Recommendations  for  Further  Work 
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Four  things  should  be  done  to  carry  this  work  forward:  1)  the 
continuum  equations  should  be  generalized  to  include  the 
convection  of  momentum  terms,  udu/dx,  to  allow  the  collector 
sheath  failure  to  be  overcome,  2)  a  collisional  transition  region 
between  the  neutral  plasma  and  the  collisionless  sheath  should  be 
developed  to  determine  the  appropriate  cut-off  for  the  plasma  ion 
distribution,  3)  the  collisional  trapping  mechanism  should  be 
explored,  and  4)  a  fast  comprehensive  sheath  model  should  be 
developed  and  inserted  into  the  implicit  non-isothermal 
computational  scheme. 
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APPENDIX  A:  ASYMTOTIC  EXPANSIONS 


A.l  Particle  Density  in  an  Accelerating  Potential 
A. 2  Particle  Density  in  an  Decelerating  Potential 


In  this  appendix  we  develop  the  asymptotic  expansions  for  the 
integrals 


These  are  needed  (  in  chap.  3  )  to  develop  the  Bohm  criterion 
correctly  and  to  demonstrate  the  necessity  of  a  general  (  as 
opposed  to  a  local  )  matching  condition.  The  expansion  of  these 
integrals  is  non-trivial  because  the  taylor  expansions  of 

^  I  U 

\fuF7T1  0n  \Tu*-X^ 

are  not  uniformly  valid  at  u  =  0  and  u  =  V%  respectively; 
therefore  we  cannot  use  the  simple  method  of  expanding  the 


-2- 


APPENDIX  A 


integrand . 


A.l  Particle  Density  in  an  Accelerating  Potential 
The  expansion  of  this  integral. 


is  done  by  splitting  the  domain  of  integration  at  a  small  fixed 
non-zero  point,  k.  On  the  lower  part  of  the  domain  the  function, 
f(u),  is  expanded  into  a  taylor  series  (  for  the  distribution 
functions  of  interest  this  is  always  possible  )  and  on  the  upper 

part  of  the  domain,  u/Au  +  X)  is  expanded  (  uniformly  because  k 
is  greater  than  zero  ) : 


The  expansion  of  the  upper  part  is  complete  and  we  need  to 
determine  the  asymptotic  expansion  of  the  lower  part.  We  will 
find  the  asymptotic  expansions  of  the  infinite  series  of 
integrands  and  regroup  the  terms  by  order  in  X  and  then  show  that 
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the  result  is  independent  of  k  as  it  should  be.  Finally,  we  will 
take  the  limit  as  k  goes  to  zero  and  we  will  then  have  the 
result. 


The  asymptotic  expansions  of  the  infinite  series  can  be  found 
recursively  starting  with  the  first  four  integrals: 

f  ufdU-  L?  ,  < 

=  T+  +  ljLr)+  0(xy 


A  \J7FfT  -T  -  + 


3  z 


OCX  ) 


*  Ofy’lnt)  . 


The  remainder  of  the  terms  are  found  by  integration  by  parts: 


,\[kv+X  ’ 


/V*-  'A. 

Cr‘-*)%v 


which  results  in 


vr 

m  25 


/W-  =JC  _  vjLl  ^>5 

^  sfu^x*  m  zcrn-z)  >  ~  * 

Adding  and  regrouping  these  terms  produces: 
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Now  we  demonstrate  that  this  result  is  independent  of  k.  The 
only  term  in  eq.  8  containing  k  is  the  X  term;  therefore  we  take 
the  derivative  of  this  term  with  respect  to  k: 


[is,  Ad- if)- 


(9) 


Since  f(u)  is  assumed  to  be  expandable  as  a  taylor  series  on 
(0,k)  the  expression  in  brackets  is  zero. 
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We  now  take  the  limit  of  eq.  8  as  k  goes  to  zero.  Through 


integration  by  parts 


therefore  we  have  the  final  result: 


lAu3mt,  = 

Jo 

+  thxfipj 

L  4  J  (id 

+  l  J hi)  uhbcJu 1 

+  0(xV. 

1/2 

The  x  term  and  the  xlnx  term  in  eq.ll  depend  only  on  the 
distribution  function  zero.  If  the  distribution  function  has  its 

1/2 

low  energy  tail  cut  off  (  as  in  chap.  3  ),  then  the  X  and  xlnx 
terms  are  zero  and  the  x  term  can  be  reduced  (  by  integration  by 

parts  )  to: 


/  **vf&  ■ 


-  X 


U,.t 


+  0  (xV . 
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A. 2  Particle  Density  in  an  Decelerating  Potential 
The  expansion  of  this  integral, 

QO 

FM  =  f  -feu)  «) 

Jjz  vul-x 

is  done  by  splitting  the  domain  of  integration  at  a  small  fixed 
non-zero  point,  k.  On  the  lower  part  of  the  domain  the  function, 
f(u),  is  expanded  into  a  taylor  series  (  for  the  distribution 
functions  of  interest  this  is  always  possible  )  and  on  the  upper 

part  of  the  domain,  u/i/(u  -  X)  is  expanded  (  uniformly  because  k 
is  greater  than  zero  ) : 


+  [  Uu)(,+  f;wt+  0 


The  expansion  of  the  upper  part  is  complete  and  we  need  to 
determine  the  asymptotic  expansion  of  the  lower  part.  We  will 
find  the  asymptotic  expansions  of  the  infinite  series  of 
integrands  and  regroup  the  terms  by  order  in  X  and  then  show  that 
the  result  is  independent  of  k  as  it  should  be.  Finally,  we  will 
take  the  limit  as  k  goes  to  zero  and  we  will  then  have  the 


result . 
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The  asymptotic  expansions  of  the  infinite  series  can  be  found 
recursively  starting  with  the  first  four  integrals: 

f  uda  __  % 

f k  uzc/u  _  ,  i  I  * 

■WVtf'-x'  T  "  -~4~ -%(i 0(f), 
CjldsL-  k.'  i  (5> 


f  V^M.  _  Li  ,, 

Jj?\nfcr  4f  +  oCx‘m). 


The  remainder  of  the  terms  are  found  by  integration  by  parts: 


-x 

<VV 


at-' 

X/*  Vv 


,  ..PT  .. 

(v*tX)**<hr  m>S 


which  results  in 


/  *3^  =  +  ±1!  x  +  a6c1A  <5> 

Y?  M  &(**) 

Adding  and  regrouping  these  terms  produces: 

J  #«><*&  -  y-r'%)  fk  u"”du. 

K  £  n!  w 


+  I  /+  +o(t‘)j 


8- 


APPENDIX  A 


+  r’w[f  -  -  X  (4- 

po  ^ 

"**  Jk^w^W’  4  ^  f  lip-  ^  *+■  0  6tV 


(7) 


Now  we  demonstrate  that  this  result  is  independent  of  k.  The 
only  term  in  eq.  8  containing  k  is  the  X  term;  therefore  we  take 
the  derivative  of  this  term  with  respect  to  k: 


=  i  r  ?r\»)L' 

JLlc^t  m  ■ 


(9) 


Since  f(u)  is  assumed  to  be  expandable  as  a  taylor  series  on 
(0,k)  the  expression  in  brackets  is  zero. 

We  now  take  the  limit  of  eq.  8  as  k  goes  to  zero.  Through 


integration  by  parts 
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-AO 

2*k 


2!*-  0 


nn.  fit 


(10) 


I  f"(u)ul»tudul 


therefore  we  have  the  final  result: 

-  fa*  «'«/-'■?■' 7 

+•  ^  +j  2  fmkiu  hj  (U) 


The  xlnx  term  in  eq.ll  depends  only  on  the  distribution  function 
zero.  If  the  distribution  function  has  its  low  energy  tail  cut 
off  (  as  in  chap.  3  ),  theT  the  xlnx  term  is  zero  and  the  x  term 
can  be  reduced  (  by  integration  by  parts  )  to: 
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APPENDIX  B:  ISOTHERMAL  PROGRAMS 

B.l  ISOTHER:  Computes  the  Isothermal  Solution 

B.2  SHIFT:  Computes  the  General  Sheath  Solution 

B.3  CAL3:  Computes  the  Sheath  Solution  given  Ushift 

and  Computed  the  Sheath  Solution  for  no 
Trapped  and  no  Surface  Emission 


-2- 


APPENDIX  B 


B.l  I SOTHER :  Computes  the  Isothermal  Solution 


ISOTHER:PROC  OPTIONS (MAIN) ; 

DCL  (T , DELT , X , DELX , FT , FX , DF IDT , DF 1DX , DF2DT , DF2DX , 

DELTIN , DELXIN , ERSH , NEWSH , DUM , 

JR ,TE , PHIW ,NCS ,NCSP ,NP0 , SURMULT, 

B 1 , DBDT , CHITEM , UZERO , SMR , LOSS ION , 

CCHI ,ACHI ,0CHI ,A1A0 ,TAU,RTAU ,TAUIN,N0 ,N1 , 

DELCHI I , DELCHIF , DELCHI S , VD , VOUT , JJR , RATIO , 

TIN , XIN , BOTEM , B ITEM , QTEM , QCTEM , RTEM , ATEM , CTEM , ABTEM , 

AO , A1 , VION , SMALL J , DELCHI .DELTAT , DELTAX , 

ERROR 1.ERROR2)  FLOAT  BIN(31); 

DCL  (FLAGCOL , I SOL , I PLASMA , I PLATE , I SINGLE , I SPEED, 

ISUR)  EXT  FIXED  BIN(31); 

DCL  (I ,J, ICASE, ICASEST, INERTIA, IDELCHI , IONIZA, IA1A0)  FIXED  BIN(31); 
DCL  (CUTOFF, TRAPPED .SURFACE ,UIN , VIONC ,UCSHIFT,CCUT, 

MU I , MUEA , KN , LAMDAR)  EXT  FLOAT  BIN(31); 

DCL  (ABINV,A,RQ,Q,BETAO ,BETA1 ,SHIFT,INVB1 )  EXT  ENTRY 
RETURNS (FLOAT  BIN(31)  ); 

CCUT=CUTOFF; 

SMR=l./4 92.2; 

IPLASMA=O;IPLATE=0; 

SURMULT® 1 . ; 

GET  FILE(ISODAT)  DATA; 

ICASE=ICASEST-1; 

LOOP : DO  DELCHI=DELCHII  TO  DELCHIF  BY  DELCHI S ; 

ICASE=ICASE+1 ; 

SHIFTLP : DO  J=1  TO  3; 

IF(FLAGCOL®0)  THEN  BEGIN; 

/*  SET  INITIAL  ITERANTS  IF  FIRST  CASE  */ 

IF  (ICASE=ICASEST)  THEN  BEGIN; 

T=TIN;X=XIN ; 

TAU=TAUIN ; 

END; 

FT=F1 (T,X) ;FX=F2(T,X) ; 

CON: DO  1=  1  TO  20; 

DELTIN=.05; 

DELXIN®. 05; 

DELT® (EXP (DELTIN )- 1. )*T; 

DELX= (EXP (DELXIN) - 1 . )*X ; 

DF1DT=(F1(T+DELT,X) -FT) /DELTIN; 

DF1DX=(F1(T,X+DELX) -FT) /DELXIN; 

DF2DT=(F2 (T+DELT ,X) -FX) /DELTIN ; 

DF2DX® ( F2 (T , X+DELX ) -FX ) /DELXIN ; 

DELTAT® (FX*DF1DX-FT*DF2DX) / (DF1DT*DF2DX-DF1DX*DF2DT) ; 

DELTAX®  ( FT*DF2DT-FX*DF1DT )  /  ( DF  1  DT*DF 2  DX  -  DF  1  DX*DF2DT ) ; 

T=EXP (DELTAT ) *T ; 


X=EXP (DELTAX)*X ; 

FT=F1(T,X); 

FX=F2(T,X) ; 

IF(ABS(FT)<=ERR0R1  &  ABS(FX)<=ERR0R2)  THEN  GO  TO  FIN; 
END  CON; 

PUT  FILE(ISOPRT)  SKIP  LIST( 'FAILED  TO  CONVERGE'); 

PUT  FILE(ISOPRT)  SKIP  DATA; 

STOP; 

END; 

ELSE  BEGIN; 

T=0.0; 

/*  SET  INITIAL  ITERANTS  IF  FIRST  CASE  */ 

IF  (ICASE=ICASEST)  THEN  BEGIN; 

X=XIN ; 

TAU-TAUIN ; 

END; 

FX=F2(T,X) ; 

CON2 :DO  1=  1  TO  20; 

DELXIN=.05; 

DELX=(EXP(DELXIN) -1 . )*X ; 

DF2DX= (F2 (T , X+DELX ) -FX ) /DELXIN ; 

DELTAX=-FX/DF2DX ; 

X=X*EXP (DELTAX ) ; 

FX=F2(T,X) ; 

IF (ABS (FX) <=ERR0R2 )  THEN  GO  TO  FIN; 

END  CON2 ; 

PUT  FILE(ISOPRT)  SKIP  LIST( 'FAILED  TO  CONVERGE'); 

PUT  FILE(ISOPRT)  SKIP  DATA; 

STOP; 

END; 

FIN: 

IF (FLAGCOL*,=0)  THEN  BEGIN; 

T=. 001 ; 

B1TEM=-ATEM*C0S (ATEM+CTEM) /SIN (ATEM+CTEM) ; 
IF(FLAGC0L-=2)  THEN  BEGIN; 

B 1 =BETA 1 (TAU , T , QCTEM ) - B ITEM ; 

C0N3 : DO  I  =  1  TO  20; 

DBDT= ( BETA 1 (TAU , T+DELT , QCTEM ) - B ITEM - B 1 ) /DELT ; 
DELTAT=-B1/DBDT; 

T=T+DELTAT ; 

B 1=BETA 1 (TAU , T , QCTEM ) - B ITEM ; 

IF(ABS(B1) <=ERR0R2 )  THEN  GO  TO  FIN1 ; 

END  C0N3 ; 

PUT  FILE(ISOPRT)  LIST(’B1  LOOP  FAILED  TO  CONVERGE'); 
STOP; 

END; 

ELSE  BEGIN; 

DUM=INVB 1 (TAU , B ITEM ) ; 

END; 

END; 

FIN1 : 

NEVSH=SHIFT (TAU , X , CHITEM .TRAPPED , SURFACE ) ; 
ERSH=NEWSH-UIN ; 


I 
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UIN , UZERO=NEWSH ; 

END  SHIFTLP;  » 

PUT  FJLE(ISOPRT)  SKIP  DATA; 

PUT  FILE (ISOPLOT)  SKIP(2)  EDIT( 'DELCHI ( ' , ICASE,  - 

')=' .DELCHI, ', ') 

(A(7),F(2),A(2),E(13,5),A(1)); 

PUT  FILE (ISOPLOT)  EDIT( *R( ICASE,  > 

•  )='  ^XEM , 1  1 ) 

(A(2),F(2),A(2),E(13,5),A(1)); 

PUT  FILE (ISOPLOT)  EDIT ( ' MACHE ( ' , ICASE , 

' )=' ,QTEM , ' , 1 ) 

(A(6),F(2),A(2),E(13,5),A(1)); 

PUT  FILE (ISOPLOT)  SKIP  EDIT( 'BETA0( ' , ICASE ,  I 

' )=' ,B0TEM, ' , 1 ) 

(A(6),F(2),A(2),E(13,5),A(1)); 

PUT  FILE(ISOPLOT)  EDIT( 'BETA1( ' .ICASE, 

')=' , B ITEM , ',') 

(A(6),F(2),A(2),E(13,5),A(1)); 

PUT  FILE (ISOPLOT)  EDIT( 'A( ICASE,  » 

' )=' ,ATEM, ' , 1 ) 

(A(2),F(2),A(2),E(13,5),A(1)); 

PUT  FILE ( ISOPLOT)  SKIP  EDIT ( ' SMALLJ ( ' , I CASE , 

' 1= 1  CM AT T 1  1  1 1 

(A(7),F(2),A(2),E(13,5),A(1)); 

PUT  FILE (ISOPLOT)  EDIT (' ECHI (', ICASE ,  | 

' )=’ ,X, ’ , ' ) 

(A(5) ,F(2) ,A(2) ,E( 13 ,5) ,A(1)) ; 

PUT  FILE (ISOPLOT)  EDIT( 'TAU( ICASE , 

(A(4),F(2)!a(2),E(13,5),A(D); 

VD=-2*(TAU-1 . ) /SMALLJ;  i 

IF(IONIZA-=0)  THEN  VD=VD-LOSSION ; 

PUT  FILE(ISOPLOT)  SKIP  EDIT( 'VD( ICASE , 

' )=' ,VD, ' , ' ) 

(A(3),F(2),A(2),E(13,5),A(1)); 

JJR=SMALLJ*EXP(-X-DELCHI ) ; 

PUT  FILE (ISOPLOT)  EDIT( 'JJR( ICASE ,  » 

• )=  *  JJR  '  '  ) 

(A(4),F(2)!a(2),E(13,5),A(1)); 

VOUT=X+DELCH I +VD ; 

IF(T<0.)  THEN  VOUT=VOUT+T ; 

PUT  FILE (ISOPLOT)  EDIT( 'VOUT( ICASE, 

' )=' ,VOUT , ' , ' )  g 

(A(5),F(2),A(2),E(13,5),A(1)); 

ACH I =TAU*LOG ( S I N ( ATEM+CTEM ) / S I N ( CTEM ) ) ; 

PUT  FILE (ISOPLOT)  SKIP  EDIT( 'ACHI (', ICASE , 

' )=' ,ACHI , '  ,  '  ) 

(A(5),F(2),A(2),E(13,5),A(1)); 

OCHI=SMALLJ*RTEM;  _ 

PUT  FILE (ISOPLOT)  EDIT( ’ OCHI (’, ICASE ,  ’ 

' )=' ,OCHI ,  ' , ' ) 

(A(5),F(2),A(2),E(13,5),A(1)); 
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PUT  FILE (ISOPLOT)  EDIT( *CCHI ( ' , ICASE , 

' )=' txt ' > ) 

(A(5),F(2),A(2),E(13,5),A(1)); 

PUT  FILE (ISOPLOT)  SKIP  EDIT( ’AO( ’ , ICASE , 

,)-,,Ao,V) 

(A(3),F(2),A(2),E(13,5),A(1)); 

PUT  FILE (ISOPLOT)  EDIT( *A1 (', ICASE, 

' )=' ,A1 , ' , ' ) 

(A(3),F(2),A(2),E(13,5),A(1)); 

PUT  FILE(ISOPLOT)  EDIT( 'C( ’ .ICASE, 

')=' , CTEM , ' , ' ) 

(A(2),F(2),A(2),E(13,5),A(1)); 

PUT  FILE (ISOPLOT)  SKIP  EDIT( 'VAVE( ' , ICASE , 

' )=' ,VION, ' , ' ) 

(A(5),F(2),A(2),E(13,5),A(1)); 

PUT  FILE(ISOPLOT)  EDIT(’UZERO( ’ .ICASE, 

')=’ .UZERO, ' , ' ) 

(A(6),F(2),A(2),E(13,5),A(1)); 

PUT  FILE(ISOPLOT)  ED IT ( 'CUTOFF ( ’.ICASE, 

1 )= 1 .CUTOFF , ' , ' ) 
(A(7),F(2),A(2),E(13,5),A(1)); 

PUT  FILE (ISOPLOT)  SKIP  EDIT ( ' TRAPPED ( ' , I CASE , 
’)=' .TRAPPED, ' , ' ) 
(A(8),F(2),A(2),E(13,5),A(1)); 

PUT  FILE (ISOPLOT)  EDIT ( ' SURFACE ( ICASE , 

’)=' .SURFACE, ' , ') 

(A(8) ,F(2) ,A(2) ,E ( 13 , 5 ) , A ( 1 ) ) ; 

PUT  FILE(ISOPLOT)  EDIT( ' 1S0L( ICASE , 

1 )=' , ISOL  ' , ' ) 

(A(5),F(2),A(2),F(2).A(1)); 

PUT  FILE (ISOPLOT)  EDIT( ' IONIZA( ICASE , 

' 1='  Tn\T7A  '  ' 1 
(A(7),F(2),A(2),F(2),A(1)); 

N0=1 . /QTEM; 

PUT  FILE(ISOPLOT)  SKIP  EDIT( ' NO (', ICASE , 

')=' .NO,*,’) 

(A (7) ,F(2) , A ( 2 ) ,E(13,5) ,A(1)) ; 

PUT  FILE (ISOPLOT)  EDIT ( ’VIONC( ICASE , 

' )= ' , VIONC , '  ' ) 

(A(6),F(2)*A(2),E(13,5),A(1)); 

PUT  FILE (ISOPLOT)  EDIT( ’UCSHIFT( ’, ICASE , 

’)=' .UCSHIFT, ' , ’) 

(A(8) ,F(2) ,A(2) ,E ( 13 , 5 ) , A ( 1 ) ) ; 

Nl=l . /QCTEM ; 

PUT  FILE (ISOPLOT)  SKIP  EDIT( ' N1 (', ICASE , 

’ )=’ ,N1 , ’ , ’ ) 

(A (7 ) ,F(2) ,A(2) ,E(13,5)  , A ( 1 ) )  ; 

PUT  FILE (ISOPLOT)  EDIT ( 'ERSH( ', ICASE , 

1 )= 1  ERSH  *  1 ) 

(A(5),F(2),A(2),E(13,5),A(1)); 

PUT  FILE (ISOPLOT)  SKIP  EDIT( ' FLAGCOL( ICASE , 
' )=' .FLAGCOL, ','  ) 
(A(8),F(2),A(2),E(13,5),A(1)); 

END  LOOP: 


RETURN; 

F1:PR0C(CX,XX)  RETURNS (FLOAT  BIN(31)); 

DCL  (CX,XX,F1X)  FLOAT  BIN(31); 

IF(IDELCHI=0)  THEN  CHITEM=0.; 

ELSE  CHITEM=DELCHI; 

QTEM=Q (TAU , XX , CH1TEM , SMALLJ , VI ON , UZERO ) ; 

/*  SURFACE  EMISSION  */ 

TE=1500. ;PHIW=1 .60; JR=20. ; 

J  JR=SMALL J*EXP (-XX-CHITEM) ; 

NCS=6 . 44E+15* ( 1500 . /TE ) ; 

NCSP=NCS/2 . / ( 1 .+2 .*EXP ( (3 . 89-PHIW)*l 1600 . /TE) ) ; 

NP0=  ( JJR*JR)*2 . 925E+1 1/QTEM* ( 1500 . /TE) ; 
SURFACE=NCSP/NP0*2 . *EXP (CHITEM ) ; 

SURF  ACE=SURF  ACE*SURMULT ; 

IF  ( I SUR=0 )  THEN  SURFACED.; 

B0TEM=BETA0 (TAU , XX , QTEM , VI  ON ) ; 

QCTEM=SQRT(2./3. 14159 )*EXP( -CX/TAU) /Cl (CX/TAU) ; 

IF (FLAGCOL=0 )  THEN  BEGIN; 

B 1TEM=BETA1 (TAU , CX , QCTEM ) ; 

END; 

ATEM=A ( BOTEM , B ITEM , QTEM , CTEM ) ; 

TAU=AB I NV ( ATEM ) ; 

A0=SQRT(3. 1415927/2 . )*QTEM*( I . /SMALLJ- 1 . )*EXP (XX/TAU) ; 
IF (CX/TAU<=0 . )  THEN 
Al=l .  ; 

ELSE 

Al=l. /Cl (CX/TAU); 

A1A0=A1/A0; 

IF  (IA1A0=0)  THEN  A1A0=1 . ; 

F 1X=CX-XX-TAU*L0G ( S I N ( ATEM+CTEM ) / S I N ( CTEM ) ) 

-TAU*LOG  ( A1A0 )  -TAU'-'LOG  ( 1 .  /  SMALLJ -1 .  )  ; 

RETURN (FIX) ; 

END  FI; 

F2 : PROC (CX ,XX)  RETURNS (FLOAT  BIN(31)); 

DCL  (CX,XX,RQTEM,F2X)  FLOAT  B IN ( 31 ) ; 

IF (  IDELCHI=0)  THEN  CHITEM=0.; 

ELSE  CHITEM=DELCHI ; 

QTEM=Q (TAU , XX , CHITEM , SMALLJ , VI ON , UZERO) ; 

/*  SURFACE  EMISSION  */ 

TE=1500 . ;PHIW=1 . 60 ; JR=20 . ; 

JJR=SMALLJ*EXP(-XX-CHITEM); 

NCS=6 . 44E+15*( 1500 . /TE) ; 

NCSP=NCS/2  .  /  ( 1 .  +2 . ’-EXP (  (3 . 89-PHIV)*l  1600  .  /  £ )  )  ; 

NPO=  ( JJR*JR)*2 . 925E+1 1/QTEM* ( 1500 . /TE ) ; 
SURFACE=NCSP/NP0*2.*EXP (CHITEM) ; 

SURF  ACE=SURF  ACE*SURMULT ; 

IF  ( I SUR=0 )  THEN  SURFACED.; 

BOTEM=BETAO (TAU , XX , QTEM , VION) ; 

QCTEM=SQRT ( 2 . / 3 . 14 159 )*EXP( -CX/TAU) /Cl (CX/TAU) ; 
IF(FLAGCOL=0 )  THEN  BEGIN; 
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B 1TEM=BETA 1 (TAU , CX , QCTEM ) ; 

END; 

ATEM=A(B0TEM,B1TEM,QTEM,CTFM); 

TAU=ABINV(ATEM) ; 

A0=SQRT(3. 1415927/2. )*QTEM*( 1 . /SMALLJ-1 . )*EXP(XX/TAU) ; 
RQTEM=RQ (ATEM , CTEM) ; 

RTEM=RAT I 0*RQTEM*QTE M / SMALL J*TAU* 

4 . /SQRT(2 .*3 . 1415927) ; 

IF(CX/TAU<=0. )  THEN 
Al=l .  ; 

ELSE 

Al=l . /C1(CX/TAU) ; 

A1A0-A1/A0 ; 

IF(IA1A0=0)  THEN  A1A0=1.; 

F2X=  2.*(1 . -TAU)/ SMALL J  +  SMALLJ*RTEM 

+TAU*LOG ( S IN ( ATEM+CTEM ) / S IN ( CTEM ) ) +XX -CX ; 

IF ( INERTI A**=0 )  THEN  BEGIN; 

F2X=F2X 

-TAU*QTEM**2/2 .*(1 . - (SIN (CTEM) /SIN (ATEM+CTEM ) )**2) ; 
END; 

IF(IONIZA-=0)  THEN  BEGIN; 

LOS  S I ON-SMR* ( V I ON+ V I ONC*S I N ( ATEM+CTEM ) / S I N ( CTEM ) ) / 
QTEM*3. 896/8. 609E-05/TE; 

F2X=F2X-LOSSION ; 

END; 

RETURN (F2X) ; 

END  F2 ; 

Cl : PROC (XX)  RETURNS (FLOAT  BIN(31)); 

DCL  (XX,C1X)  FLOAT  BIN(31); 

IF (XX<=0 . )  THEN  C1X=1 . ; 

ELSE  C1X=1.+ERF(SQRT(XX)); 

RETURN (C IX); 

END  Cl; 

C2 : PROC (XX)  RETURNS (FLOAT  BIN(31) ) ; 

DCL  (XX,C2X)  FLOAT  BIN(31); 

C2X=EXP(XX)*(1 . -ERF(SQRT(XX) ) ) ; 

RETURN (C2X); 

END  C2; 

END  ISOTHER; 
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B.2  SHIFT :  Computes  the  General  Sheath  Solution 


SHIFT : PROC (CE , PHIM , R .TRAPPED , SURFACE )  RETURNS (FLOAT  BIN(31)) 
DCL(BETA,BETATEM,CI,SCI,CE,R,UR, 

NEB.NPLASMA, 

Cl ,C2,C5 ,C6 ,PHIM, 

PLASMA, SURFACE, 

CHI, 

URL,URH,MURM,URM, 

TEM1 ,TEM2,TEM3 ,TEM4 ,TEM5 ,TEM6 ,TEM7) 

FLOAT  BIN(31) ; 

DCL(  I , ICASE)  FIXED  BIN(31); 

DCL(FIB , FID, FI I ) 

EXT  ENTRY (FLOAT  BIN(31) .FLOAT  BIN(31), 

FLOAT  BIN(31), FLOAT  BIN(31) .FLOAT  BIN(31)) 

RETURNS (  FLOAT  BIN(31)); 

DCL  (FEB.FIJ)  EXT  ENTRY (FLOAT  BIN(31) .FLOAT  BIN(31), 
FLOAT  BIN(31), FLOAT  BIN(31)) 

RETURNS (  FLOAT  BIN(31)  ); 

DCL  (TRAPPED)  FLOAT  BIN(31); 

DCL  (M)  FIXED  BIN(31)  EXT; 

M=20 ; 

CI=1 . ;BETATEM=. 1 ; 

SCI=SQRT(CI); 

BETA=BETATEM*SQRT(CE/2 . ) ; 

URH=3. ;URL=-3. ; 

LOOPUR : DO  I  =  1  TO  10; 

URM= (URL+URH ) / 2 . ; 

MURM=M I N I NT ( URM ) ; 

IF (MURM<=-0. 00000001)  THEN  URL=URM ; 

ELSE  URH=URM ; 

END  LOOPUR; 

FINUR : UR=URM*SQRT(2 . /CE) ; 

RETURN (UR); 

MININT : PROC (URTEM)  RETURNS (  FLOAT  BIN(31 ) ) ; 

DCL  ( URTEM, FNETM IN, TEMP)  FLOAT  BIN(31); 

DCL  J  FIXED  BIN(31) ; 

CHI=PHIM;NEB=(EA(URTEM,R) -K3 (CHI ) ) /K4 (CHI ) ; 

CHI=0 . ;FNETMIN=EA (URTEM ,R) -K3 (CHI ) -K4 (CHI )*NEB ; 
LOOP: DO  J=1  TO  25; 

CHl=J*PHIM/25 . ; 

TEMP=EA(URTEM,R)-K3(CHI)-K4(CHI)*NEB; 

FNETM I N=M I N ( FNETM I N , TEMP ) ; 

END  LOOP; 

LOOP 2 : DO  J=1  TO  10; 

CHI=J*PHIM/250. ; 

TEMP=EA (URTEM ,R) -K3 (CHI ) -K4 (CHI )*NEB ; 

FNETM I N=M I N ( FNETM I N , TEMP ) ; 

END  LOOP 2 ; 
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RETURN (FNETMIN); 

END  MININT; 

EA : PROC (TUZERO , TR )  RETURNS (FLOAT  BIN(31)); 

()CL  (TUZERO ,  TR ,  TEMEA ,  TPHIM ,  TPHI Z ,  SQZM ,  SQPI , 

IOXZM , IZMXX , IOM , 

SQZX.IZXI ,  IOX, IZMXI ,SQXZM,IZI ,IZMI) 

FLOAT  BIN(31) ; 

TPHIM=PHIM;TPHIZ=(PHIM+TR) ; 

IF (TPHIZ>=TPHIM)  THEN  BEGIN; 

SQZM=SQRT(TPHIZ-TPHIM) ; 

SQZX=SQRT (TPHI Z -CHI ) ; 

IZMI=(1 . -ERF (SQZM) )*SQRT(3 . 1415926)/2 . ; 
IOM=ERF(SQRT(TPHIM) )*SQRT(3. 1415926)/2. ; 

IZXI=( 1 . -ERF(SQZX) )*SQRT(3 . 14 15926 ) /2 . ; 

I OX=ERF ( SQRT (CHI ) )*SQRT(3 . 1415926 )/2 . ; 

IZI=(1 . -ERF (SQRT (TPHIM) ) )*SQRT(3 . 1415926 )/2 . ; 

IZMXI=(1 . -ERF (SQRT (TPHI Z-TPHIM+CHI ) ) )*SQRT(3 . 1415926 ) /2 . 
PLASMA= ( SQRT ( 3 . 1 4 1 5  9  2  6 ) - SURFACE* I ZMI ) / 

(FID(BETA/SCI, 1000000. ,1. .TUZERO, 0. )  + 

FID (BETA/SCI, SQZM, 1. .TUZERO, 0.)  ); 

TEMEA=2 . *PLASMA* 

(  FI I (BETA/SC 1, 1000000. ,1. .TUZERO, 

SQRT(CHI))  - 

FII (BETA/SCI , 1000000 1 . .TUZERO , 

0.)  + 

FII (BETA/SCI, SQZM, 1. .TUZERO, 

SQRT(CHI))  - 

FII (BETA/SCI , SQZM, 1 . .TUZERO, 

0.)) 

+  TRAPPED* 

(  2.*EXP(CHI)*I0X-2.*SQRT(CHI)  ) 

+SURFACE* 

(  EXP (TPHIM -TPH I Z)*SQRT (TPHI Z-TPHIM+CHI ) 

+EXP(CHI )*IZMXI 

-SQRT (TPHI Z -TPHIM )*EXP(TPHIM-TPHIZ) 

-IZMI  ); 

TEMEA=TEMEA/ SQRT ( 3 . 1415926) ; 

END; 

ELSE  BEGIN; 

SQPI=SQRT(3. 1415926); 

PLASMA=SQPI*(1. -SURFACE/2. )/ 

FID(BETA/SCI, 1000000. , 1 . .TUZERO, 0 . ) ; 

TEMEA= 

2. *PLASMA*(FI I (BETA/SCI, 1000000. ,1. .TUZERO, 
SQRT(CHI))- 

FI I (BETA/SCI, 1000000. ,1. .TUZERO, 

0.)  )/SQPI; 

IF (CHI<=TPHIM-TPHIZ)  THEN  BEGIN; 

I 0X=ERF ( SQRT (CHI))*SQPI/2. ; 

TEMEA=TE.MEA  +  SURFACE* 
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(SQPI/2.*(EXP(CHI)-1. ) 
+EXP(CHI)*IOX-SQRT(CHI))/SQPI ; 

END; 

ELSE  BEGIN; 

SQXZM=SQRT(CHI+TPHIZ-TPHIM) ; 

I0XZM=SQPI*ERF(SQXZM)/2. ; 

I ZMXI=( 1 . -ERF ( SQRT (TPHI Z -TPHIM+CHI ) ) )*SQPI / 2 . ; 
IZMXX=(  ERF(SQRT(CHI)) 

-ERF (SQRT (TPHI Z -TPHIM+CHI ) )  )*SQPI/2. ; 
TEMEA=TEMEA  +  SURFACE* 

(SQRT(TPHIZ-TPHIM+CHI )*EXP(TPHIM-TPHIZ)*2 . 
+EXP (CHI )*IZMXI 

-SQPI/2. -SQRT(CHI)+EXP(CHI)*IZMXX)/SQPI 
+TRAPPED* (2 .*EXP (CHI )*IOXZM  - 

2.*EXP(TPHIM-TPHIZ)*SQXZM)/SQPI; 


END; 


END; 

RETURN (TEMEA*CI ) ; 

END  EA; 

K3 : PROC (CHITEM)  RETURNS (  FLOAT  BIN(31) ) ; 

DCL(CHITEM)  FLOAT  BIN(31); 

TEMl=SQRT(3.1415926)/2. ; 

TEM3=SQRT(3 . 1415926)/2 .*ERF(SQRT(PHIM/CE) ) ; 
C1=(TEM1+TEM3)/ (TEM1) ; 

C5=  (l.-EXP(-CHITEM/CE))+  EXP(-CHITEM/d.  )*( 

2*FI 1(0. ,SQRT((PHIM-CHITEM)/CE),1. ,0. .SQRTCCHITEM/CE) ) 
2*FI I (0 . ,SQRT((PHIM-CHITEM)/CE),1. ,0. ,0. )  )* 

2. /SQRT(3. 1415926)+ 

2*FI 1(0. ,SQRT(CHITEM/CE) , 1 . ,0. , 0 . )*2 . /SQRT(3 . 1415926) ; 

RETURN (C5*CE/C1 ) ; 

END  K3 ; 

K4: PROC (CHITEM)  RETURNS (  FLOAT  BIN(31) ) ; 

DCL(CHITEM)  FLOAT  BIN(31); 

TEM1=SQRT(3 . 1415926)/2. ; 

TEM2=TEM1 ; 

TEM3=SQRT(3. 1415926)/2 .*ERF(SQRT(PHIM/CE) ) ; 

TEM4=FID(0. ,1000000. ,1. , 0 . .SQRT(PHIM) ) ; 

TEM6=2*FII(0. ,1000000. ,1. ,0. .SQRT(PHIM) ) ; 

TEM7=2*FII(0. ,1000000. , 1 . ,0. , SQRT(PHIM-CHITEM) ) ; 


C1=(TEM1+TEM3)/(TEM1); 

C2=TEM4/TEM2 ; 

C5=  (1 . -EXP(-CHITEM/CE) )+  EXP(-CHITEM/CE)*( 

2*FII (0. ,SQRT((PHIM-CHITEM)/CE),1. ,0. ,SQRT(CHITEM/CE) ) 
2*FI 1(0. ,SQRT((PHIM-CHITEM)/CE) ,1. ,0. ,0. )  )* 

2. /SQRT(3. 1415926)+ 

2*FII(0. ,SQRT(CHITEM/CE) , 1 . ,0. ,0. )*2 . /SQRT(3 . 1415926) ; 
C6= (TEM6 -TEM7 ) /TEM2 ; 
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B.3  CAL3 :  Computes  the  Sheath  Solution  given  Ushift 

and  Computed  the  Sheath  Solution  for  no 
Trapped  and  no  Surface  Emission 


CAL3 : PROC (ECH I , PHI Z , CE , C I , BETATEM , MR , ERROR , UZERO , BNEB , BJNET , VAVE ) ; 
DCL (BETA , BETATEM , C I , SCI , CE , ALPH , BET , R , UR , 

NEB, NEP, BNEB, BJNET.JNET, VAVE, 
C1,C2,C3,C4,C5,C6,H,PHIM,ECHI, 

F0,F1 ,F2 ,UR0 ,UR1 ,UR2 , 

K1 ,K2,K3,K4, 

JION, JEE.RZ, 


JEP, PLASMA, 

DELX , DELTAX , X , HX , DHDX , 

UP,MR,ATEM, 

ERROR, PHIZ, TEM1, 

TEM2 , TEM3 , TEM4 , TEM5 , TEM6 , TEM7 , UZERO ) 

FLOAT  BIN(31) ; 

DCL (FIB , FID, FI I ) 

EXT  ENTRY (FLOAT  BIN(31) .FLOAT  BIN(31), 

FLOAT  BIN(31), FLOAT  BIN(31) .FLOAT  BIN(31)) 

RETURNS (  FLOAT  BIN(31)); 

DCL  (FEB.FIJ)  EXT  ENTRY (FLOAT  BIN(31) .FLOAT  BIN(31), 

FLOAT  BIN(31), FLOAT  BIN(31)) 

RETURNS (  FLOAT  BIN(31)  ) ; 

DCL(I , J)  FIXED  BIN(31) ; 

DCL (M, I SOL, I  SINGLE , I SPEED , I PLASMA , IPLATE)  FIXED  BIN(31)  EXT; 
DCL  (TRAPPED, SURFACE, INTERR,NENO,UIN)  FLOAT  BIN (31)  EXT; 


j  *************^V*-^****--'n’rVr******-!Wr**-.'r******-,V**T,r********  j 

/*  */ 
/*  EXTERNAL  VARIABLE  INPUTS  */ 

/*  ISOL  =  0  FOR  CUTOFF  ION  DISTRIBUTION  BOHM  */ 

/*  =  1  SQRT(CHI)  MATCHING  */ 

/*  IS INGLE =  0  FOR  DOUBLE  SHEATH  */ 

/*  =  1  SINGLE  SHEATH  */ 

/*  I SPEED  =  0  MATCHED  SHEATH  */ 

/*  =  1  VSHIFT  INPUT  */ 

/*  =2  VSHIFT  =  VAVERAGE  */ 

/*  IPLASMA=  0  FOR  FULL  PLASMA  ELECTRON  DYNAMICS  */ 

/*  IPLATE  =  0  FOR  FULL  PLATE  ELECTRON  DYNAMICS  */ 


/*  */ 
/1rMt**-Mrh-Mr^-Mit«****it*******Mt*********1li*****rMeii«****  j 

PHIM=ECHI ; 

R=(PHIZ-PHIM)/CI ; 

SCI=SQRT(CI ) ; 

BETA=BETATEM*SQRT ( CE / 2 . ) ; 


TEM1=SQRT(3. 1415926)/2.  ; 


TEM2=TEM1 ; 

TEM3=SQRT(3. 1415926)/2.*ERF(SQRT(PHIM/CE)) ; 

TEM4=FID(0. ,1000000. ,1 . ,0 . .SQRT(PHIM) ) ; 

/*  TEM4=SQRT ( 3 . 1 4 1 5 9 26 ) / 2 . *EXP ( PH 1M )* ( 1 . -ERF ( SQRT ( PHI M ) ) ) ;  * / 
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TEM5=FEB(0. ,1000000. ,1. ,PHIM); 

/*  TEM5=1 . /SQRT(PHIM) -2 .*EXP(PHIM)*SQRT(3 . 1415926)/2 .* 

(1 . -ERF(SQRT(PHIM) ) ) ; 

TEM6=2*FII(0. ,1000000. ,1. ,0. .SQRT(PHIM) ) ; 

TEM7=2*FII (0. , 1000000. , 1 ., 0.  , 0.  ) ; 

C1=(TEM1+TEM3)/(TEM1); 

C2=TEM4/TEM2 ; 

C3=2 . * (TEM1+TEM3 ) /TEM1+SQRT (CE/PHIM)*EXP ( -PHIM/CE ) /TEM1 ; 
C4=TEM5/TEM2 ; 

C5=(  TEM3 -SQRT (PHIM/CE )*EXP ( -PHIM/CE )+ 

SQRT(3 . 1415926 )/2 .*( 1 . -EXP ( -PHIM/CE))  )/ 

(TEM1) ; 

C6=(TEM6-TEM7)/TEM2 ; 

IF  (IPLASMA->=0)  THEN  BEGIN; 

Cl=2. ; 

C3=4.  ; 

C5=2 . * ( 1 . -EXP ( -PHIM/CE ) ) ; 

END; 

IF(IPLATE->=0)  THEN  BEGIN; 

C2=1./SQRT(1.+3.14159*PHIM); 

C4=3. 14159/2. / (1. +3. 14159*PHIM)**( 1.5); 

C6=2 . /3 . 14159*(SQRT(1 .+3 . 14159*PHIM) -1 . ) ; 

END; 

K1=C3/CE/C1 ; 

K2=C4+C2*C3/CE/C1 ; 

K3=C5*CE/C1 ; 

K4=C6-C2*C5*CE/C1 ; 


IF (ISPEED=0)  THEN  BEGIN; 

IF(ISOL=0)  THEN  BEGIN; 

UR 0=2 . 0 ; F0=F (URO ) ; 

UR1=2 . 1 ;F1=F(UR1) ; 

END; 

ELSE  BEGIN; 

UR0=1 . 0 ;FO=F (URO) ; 

UR1=1 . 1 ;F1=F(UR1) ; 

END; 

CON : DO  1=  1  TO  20; 

UR2=UR1- (UR1 -URO)/ (F1-F0)*F1 ; 

F2=F (UR2) ; 

IF(IS0L=0)  THEN  BEGIN; 

IF(ABS (F2 )<=ERROR)  THEN  GO  TO  FIN; 

END; 

ELSE  BEGIN; 

IF(ABS(F2) <=ERROR '^TRAPPED )  THEN  GO  TO  FIN; 
END; 

UR0=UR1 ;F0=F1 ; 

UR1=UR2 ;F1=F2 ; 

END  CON; 

PUT  FILE(AUXPRT)  SKIP  LIST( 'FAILED  TO  CONVERGE’); 
STOP; 

F : PROC(XX)  RETURNS (FLOAT  BIN(31) ) ; 
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DCL  (FX.XX)  FLOAT  BIN(31); 

IF(ISOL=0)  THEN  BEGIN; 
ALPH=(K1-BA(XX,R))/K2; 

BET= (EA(XX,R)-K3)/K4; 

IF(  ISJNGLE=0)  THEN 
FX=ALPH-BET ; 

ELSE 

FX=ALPH-NENO; 

END; 

ELSE  BEGIN; 

BETA=0 . ; 

IF(R>=0.)  THEN  BEGIN; 

FX=2 .*TRAPPED-EXP(-XX**2)/ 

(FID(BETA/SCI, 1000000. ,1. ,XX,0.)  + 
FID(BETA/SCI ,SQRT(R) ,1. ,XX,0. )  ); 

END; 

ELSE  BEGIN; 

FX=SURFACE -EXP ( -XX**2 ) / 

(FID (BETA/SCI, 1000000. ,1. , XX, 0. )  ); 

END; 

END; 

RETURN (FX) ; 

END  F; 

END; 

IF(ISPEED=1)  THEN  BEGIN; 

UR2=UIN*SQRT(CE/2 . /Cl ) ; 

END; 

IF ( ISPEED=2)  THEN  BEGIN; 

IF(R<=0.)  THEN  UR2=UIN; 

ELSE  UR2=UIN*EXP(-R)*SQRT(CE/2./CI) ; 

END; 


FIN:UR=UR2; 

IF(IS0L-=0)  THEN 

ALPH , BET= (EA (UR , R) -K3 )/K4 ; 

NEB=(BET+ALPH)/2 . ; 

IF(ISINGLE-=0)  THEN 
NEB=NEN0 ; 

INTERR=NEB-BET; 

NEP=(1 . -C2*NEB)/C1 ; 

IF  R>0 .  THEN 

UP=MAX (BETA , SQRT (R)  •'•SCI ,  0 .  )  ; 

ELSE  UP  =  MAX(BETA,0. ) ; 

IF  R<=0 .  THEN  RZ=0 . ; 

ELSE  RZ=R; 

ATEM= (  FI J(UP/SCI, 1000000. ,1. , UR)  -  SURFACE*. 5*EXP(-R)  )/ 

(  FID(BETA/SCI , 1000000 . , 1 , ,UR, 0. )  +  FID(BETA/SCI .UP/SCI , 1. ,UR,0. ) 
+  SURFACE*SQRT(3 . 1415926) /2 .*(1 . -ERF ( SQRT (RZ) ) )  ); 

JION=SQRT ( C I /MR ) *ATEM ; 

VAVE=SQRT ( 2 . *CI/CE)*ATEM ; 
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JEE=  -NEB/SQRT(3. 1415926) ; 

JEP=  - (1-C2*NEB)/C1/SQRT(3. 1415926)*SQRT(CE)* 
EXP(-PHIM/CE) ; 

JNET=- JEE+JEP; 

BNEB=N$B; 

B JNET= JNET / SQRT (CE / 2 . ) ; 

UZER0=UR*SQRT(2.*CI/CE) ; 

PUT  FILE(REFPRT)  SKIP  DATA(C1,C2,C3,C4,C5,C6,NEB) ; 

RETURN; 

BA : PROC (TUZERO , TR )  RETURNS (FLOAT  BIN(31)); 

DCL(TUZERO ,TR , SQTR ,TEMBA ,  B 1 , B2 , B3 , B4 )  FLOAT  BIN( 31 ) ; 
B1=FIB(BETA/SCI, 1000000. ,1. .TUZERO, 0. ); 
B3=CI*FID(BETA/SCI ,1000000. ,1. , TUZERO, 0. ) ; 

IF (TR> (BETA/ SCI )**2 )  THEN  BEGIN; 

SQTR=SQRT(TR) ; 

B2=FIB( BETA/SCI, SQTR, 1. , TUZERO, 0. ); 

B4=CI*FID( BETA/SCI ,SQTR, 1 . , TUZERO, 0. ), 
TEMBA=(B1+B2)/(B3+B4); 

END; 

ELSE 

TEMBA=B1/B3 ; 

RETURN (TEMBA); 

END  BA; 

EA: PROC (TUZERO, TR)  RETURNS (FLOAT  BIN(31>); 

DCL (TUZERO , TR , TEMEA , TPH I M . TPH I Z .  SQZ'I ,  SQZ ,  SQ*! .  SQP I  , 
IZI , IZMI , IOM, IOZ, IZM) 

FLOAT  BIN(31); 

TPH I M=PH I M / C I ; TPH I Z=TPH I M+TR ; 

IF (TPHIZ<=0 . )  THEN  TPHIZ=0.; 

IF(TPHIZ>=TPHIM)  THEN  BEGIN; 

SQZM=SQRT ( TPH I Z -TPH I M ) ; 

IZMI=( 1 . -ERF(SQZM) )*SQRT(3 . 14]5926)/2.  ; 
IOM=ERF(SQRT(TPHIM) )*SQRT ( 3 . 1415926  )/2 . ; 

IZI=(1 . -ERF ( SQRT (TPHIM ) ) )*SQRT(3 . 1415926  )/2 .  ; 
PLASMA= ( SQRT ( 3 . 1415926 ) -SURFACE* I ZM I )/ 

(FID(BETA/SCI, 1000000. , 1 . , TUZERO, 0. I  ♦ 
FID(BETA/SCI ,SQZM,1. , TUZERO, 0. )  ) ; 

TEMEA=2 . *PLASMA* 

(  F1I (BETA/SCI , 1000000. , 1 . .TUZERO, 

SQRT(TPHIM) )  - 

FII (BETA/SCI, 1000000. ,1. .TUZERO, 

0.)  + 

FII (BETA/SCI, SQZM.l. .TUZERO, 

SQRT(TPHIM) )  - 

FII (BETA/SCI ,SQZM,1. .TUZERO, 

0.)) 

+  TRAPPED* 

(  2.*EXP(TPHIM)*I0M-2.*SQRT(TPHIM)  ) 

+  SURFACE* 

(  EXP(TPHIM)*IZI+EXP(TPHIM-TPHIZ)*SQRT(TPHIZ) 

- I ZMI - SQZM*EXP ( TPH I M -TPH I Z )  ) ; 
TEMEA=TEMEA/SQRT(3. 1415926); 


-16-  APPENDIX  B 

END; 

ELSE  BEGIN; 

SQZ=SQRT(TPHIZ); 

SQM=SQRT(TPHIM) ; 

SQPI=SQRT(3. 1415926); 

I0Z=SQPI*ERF(SQZ)/2. ; 

IZI=SQPI/2. -IOZ; 

IZM-SQPI/2 . * (ERF (SQM) -ERF(SQZ)) ; 

PLASMA=SQPI*(1 . -SURFACE/ 2. )/ 

FIDCBETA/SCI, lOOOOOO. ,1. .TUZERO, 0.); 
TEMEA=SURFACE* (EXP (TPHIM)* ( IZI+IZM) 

+2.*SQZ*EXP(TPHIM-TPHIZ) 
-SQM-SQPI/2. )/SQPI 
+TRAPPED*(2.*EXP(TPHIM)*I0Z  - 

2 . *EXP (TPHIM-TPHIZ)*SQZ) /SQPI 
+2. *PLASMA*(FII (BETA/SCI ,1000000. ,1. .TUZERO, 
SQRT(TPHIM))- 

FII (BETA/SCI, 1000000. ,1. .TUZERO, 
0.)  )/SQPI ; 

END; 

RETURN (TEMEA*CI) ; 

END  EA; 

END  CAL3 ; 
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APPENDIX  C:  NON- ISOTHERMAL  PROGRAMS 

C.l  PREDCOR:  Computes  the  Unsteady 

Non-Isothermal  Solution 

C.2  SHEATH:  Computes  the  Sheath  Model  Solution 
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C.l  PREDCOR :  Computes  the  Unsteady  Non -Ios thermal  Solution 


/★  PREDCOR  FOR  FULL  TEC  */ 

/*  APRIL  27,  1980:  REVISED  TO  INCLUDE  YEN  FACTOR  */ 

/*  MAY  4,  1980:  REVISED  TO  INCLUDE  IVD,  FINAL  DOT,  ELOSSB,  */ 

/*  AND  CORRECT  MURS  IN  BC.  */ 

PREDCOR : PROC (  Tl,  T2,  TAU,  NEB,  NSTEPS,  TE,  TC,  ENR,  CNR, 

TDOT1,  ND0T1, 

PN,  SMR,  LAMDAR,  KN,  NR,  ARECN,  TCHAR,  EGNDB ,  RE, 

N,I,  AN,  AT,  BN,  BT,  CN,  CT)  REORDER; 

DECLARE 


SUMV  ENTRY ( (★)  FLOAT  DEC(16) .FIXED  BIN(31) .FIXED  BIN(31)) 


RETURNS (FLOAT  DEC (16)), 

SQRT  BUILTIN, 

/★THE  FOLLOWING  TIME  VARIABLES  ARE  NOND  BY  TCHAR.  */ 

(T1,T2,  /★ START  &  FINISH  TIMES.  ★/ 

DT,  /*ACTUAL  TIME  STEP  USED.  */ 

TIME, 

AN , AT , BN , BT , CN , CT ,  /*PRED-COR  ALPHA , BETA , GAMMA  */ 

IVD  EXT,  /*PLASMA  POWER  GAIN  */ 

(LAMTAU,LAMNEB)EXT, 

/★THE  FOLLOWING  VARIABLES  REFER  TO  THE  MOST  RECENT  TIME  */ 


TAU(*) ,  /*E-  TEMPERATURE  (NOND  BY  TE) .  ★/ 

NEB(*) ,  /*ELECTRON  DENSITY  (NOND  BY  NR)*/ 

TD0T1(*) ,  ND0T1(*),  /*PREDICTOR  STEP  TIME  DERIV.S  ★/ 

( (ENE , CNE )  INIT(0.8) ,  /*E  &  C  EMITTED  DENSITY,  NE.  */ 

(ECHI.CCHI)  INIT(3) ,  REMITTER  &  COLLECTOR  DROPS  */ 

(EALPHA , CALPHA)  INIT(0.5)  /*E  &  C  ION  SPEED  PARAMETERS  */ 


)  STATIC  EXTERNAL, 

/*THE  FOLLOWING  ARE  CONSTANT  DURING  THIS  PROGRAM  */ 


MUI(0:N+1) , 

/*ION  MOBILITY 

*/ 

ONE  INIT(l) , 

I, 

/*CURRENT  (NOND  BY  REF  DIFF  C) 

*/ 

DZ, 

/*ZETA  INCREMENT  BETWEEN  PTS 

*/ 

TCHAR, 

/*AN  ELECTRON  TRANSIT  TIME 

*/ 

ERR, 

/*ACCURACY  PARAM  FOR  SHEATH 

*/ 

TE ,TC, 

REMITTER  &  COLLECTOR  TEMPERATU*/ 

DTAUNDZ , 

PN, 

/*NEUTRAL  PRESSURE  (TORR) 

*/ 

ENR, CNR, 

/*E.  &  C.  RICHARSON  DENSITIES. 

*/ 

NR, 

/★REFERENCE  ELECTRON  DENSITY 

*/ 

EGNDB, 

/*E(0)/KT(E) .  NON-D  BINDING  E 

*/ 

ELOSSB  EXT, 

/★ENERGY  LOSS  PER  IONIZATION 

*/ 

RE, 

/★Q(E-A)  =Q0  *  E**-RE 

*/ 

SMR, 

/★SQRT  OF  ELECTRON/ION  MASS  RAT*/ 

LAMDAR, 

/★  MFP  RATIO,  =RMUR/SMR 

*/ 

RMUR, 

/★MU  RATIO,  =SMR*LAMDAR 

*/ 

KN, 

/★KNUDSEN  NUMBER 

*/ 
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nnr,  /^reference  neutral  density.  */ 

arecn,  /*coefficient  of  recombination  */ 

PI,  /*  3.14159...  */ 

ca.csaha,  /^constants  in  msource  eqn  */ 

/*THE  FOLLOWING  ARE  VECTORS  (0:N+1).  */ 

(NNB.TAUN,  /*NEUTRAL  DENSITY  &  TEMP.  */ 

/*  NON-D  BY  EMITTER  VALUES  */ 
MSOURCE,  /*E-  PRODUCTION  RATE  (NOND)  */ 

ESOURCE,  /*ENERGY  SOURCE  TERM  (NOND)  */ 

CV, 

muea,  /*e-  mobility  among  atoms.  */ 

NDOT2.TDOT2, 

TTILDA.NTILDA) 

(0:N+1)  )  FLOAT  DEC(16), 

CFIX  BIT(l)  EXT,  /*DID  C-SHEATH  REQUIRE  FIX?  */ 

EFIX  BIT(l)  EXT, 

(NSTEPS, 

N,  /*#  OF  GRID  PTS,  E  TO  C  INCL.  */ 

J, 

COUNT)  /^PRESENT  */ 

FIXED  BIN(31) ; 

DCL  IDEN  FIXED  BIN(31)  EXT; 

/*HANDLE  EXCEPTIONAL  CONDITIONS.  */ 

ON  FINISH  PUT  SKIP (5)  DATA; 


PI=3. 1415926  +  5.3589793E-8; 

DZ=ONE/(N-l); 

DT=(T2-T1) /NSTEPS; 

CFIX,EFIX=’0'B; 

ERR=lE-3 ; 

/*SET  NEUTRAL  TEMPERATURE  AND  DENSITY.  */ 

IF  TE=TC  THEN  TAUN=1; 

ELSE  DO  J=0  TO  N+l; 

TAUN(J)=1  +  (TC/TE-1)*(J-1)/(N-1); 

END; 

NNR=965.5E16*PN/TE; 

NNB=1/TAUN; 

DTAUNDZ=TAUN (N) -TAUN  ( 1 ) ; 

/*SET  TRANSPORT  PARAMETERS.  */ 

RMUR=LAMDAR*SMR ; 

MUI=SQRT(TAUN) ; 

/*SET  IONIZATION  AND  SAHA  PARAMETERS.  */ 

CA=0.41283*ARECN*TCHAR*(NR/1E14)**2  *  (TE/1500)**>4.5; 
CSAHA=LOG(  (1.4027E20*NNR/NR/NR)  *  (TE/1500)**1 .5  ); 

DO  COUNT=0  TO  NSTEPS- 1; 

TIME=Tl+COUNT*DT; 

/*PREDICTOR  STEP  */ 

CALL  DOT(NDOTl ,TDOTl ,NEB ,TAU) ; 

NTILDA=NEB  +  AN*DT*ND0T1; 

TTILDA=TAU  +  AT*DT*TD0T1; 

/^CORRECTOR  STEP  */ 

IF(  AN=0.  &  AT=0 .  )  THEN  BEGIN; 


( 
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ND0T2=0. ;TDOT2=0. ;G0  TO  NOCOR;END; 

CALL  DOT(NDOT2 ,TDOT2 ,NTILDA,TTILDA) ; 

NOCOR: 

NEB=NEB  +  DT*(  BN*NDOTl+CN*NDOT2  ); 

TAU=TAU  +  DT*(  BT*TDOTl+CT*TDOT2  ); 

END; 

/*UPDATE  TIME  DERIV.S,  IMAGE  POINTS,  AND  FIND  PLASMA  POWER  GAIN*/ 
CV(0),CV(N+1)=0; 

ESOURCE(O) ,ES0URCE(N+1)=0; 

CALL  DOT(NDOTl ,TDOTl ,NEB ,TAU) ; 

IVD=  0.5*(  ESOURCE ( 1 ) -CV ( 1 )*TDOTl ( 1 )  ) 

+  SUMV(  ESOURCE -CV*TD0T1,  2,  N-l) 

+0 . 5*(  ESOURCE (N)-CV(N)*TD0T1(N)  ); 

IVD=DZ*IVD+  2*I*(TAU(1) -TAU(N) )  -  (NEB(1)*ENE/KN)*(TAU(1)-1) ; 


DOT:  RETURNS  WITH  NEW  NEBDOT  AND  TAUDOT. 


/*  COLLECTOR  EMISSION  IS  NEGLECTED. 

/* 

/*  THERMAL  DIFFUSION  RATIO  IS  INCLUDED. 

/* 

/*  APRIL  2,  1980,  (B.C.  BASED  ON  MARCH  1,  1980  VERSION) 


PROC (NEBDOT,  TAUDOT,  NEB,  TAU)  REORDER; 

DECLARE 

CHKDOT  BIT(l)  EXT  INIT('O'B) ,/*IF  1  PRINT  DIAGNOSTIC  INFO  */ 
SHEATH  ENTRY (DEC (16) ,DEC(16) ,DEC(16) ,DEC(16) ,DEC(16) , 

DEC (16) , DEC (16) ,DEC(16) ,DEC(16) ,DEC(16)) , 

(TAUDOT (*), 

NEBDOT (*) , 

NEB(*) , 

TAU(*) , 

PC(0 :N+1) ,  /*CHARGED  PARTICLE  PRESSURE  */ 

F(15)  EXT  INIT(5 . 74E-3,  1.4E-3,  2.3,  .2,  .027,  .00574, 

.0424,  3.2,  61.893,  11.607,  15473,  27.04), 


/*SHEATH  VARIABLES  */ 

NCMIN,  /*MINIMUM  NEB (N)  TO  ALLOW  I.  */ 

(  U , GU , DELU , DGDU , DELTAU , 

^TRANSPORT  VARIABLES  */ 

FYEN  EXT  INIT(l) ,  /*YEN  THERMAL  CONDUCTIVITY  FACT*/ 

K(0 :N+1) ,  /*THERMAL  CONDUCTIVITY  */ 

MNS ,  /*SUM  OF  MUI*NEB  AT  J  &  J+l  */ 

MUIS ,MURS ,  /*SUM  OF  MU I  &  MUR  AT  J,J+1  */ 

'  /CONSERVATION  EQUATION  VARIABLES  */ 

EN_Z ,CN_Z ,  /*SPATIAL  NEB  DERIVS  */ 

ET_ETA , ET_Z , CT_ETA , CT_Z ,  /*SPATIAL  DERIVS  OF  TAU  */ 

EDETA , CDETA ,  /*ETA  XPACING  BETWEEN  GRID-PTS  */ 

EPC_Z,CPC_Z,  /*PC  GRADIENT  FROM  B.C.  */ 

MURSOLD .MNSOLD ,MUI SOLD , 

<•  NBCOA  ,TBC0A ,  NBCOB  .TBCOB ,  NBCOC  ,TBC0C , 

NBC 1A , TBC 1 A , NBC 1 B ,TBC1B ,NBC1C ,TBC1C , 

NEBA(0:N+1) ,NA(0:N+1) ,NEBB(0:N+1) ,NB(0:N+1) , 


( 
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NEBC(0:N+1) ,NC(0  N+l) ,NEBV(0:N+1) ,ND(0:N+1) , 

NS (0 : N+l ) , TAUA (0 : N+l ) , TA (0 : N+l ) , TAUB (0 : N+l ) , 

TB(0:N+1) ,TAUC(0 :N+1) ,TC(0:N+1) ,TAUV(0:N+1) ,TS(0:N+1) , 

NEBU(0:N+1) ,TAUU(0:N+1) , 

ALPHA (0 : N+l ) , BETA (0 : N+l ) , 

/*TEMP0RARY  ENERGY  EQUATION  VARIABLES  */ 

QKM,QKP, 

DETA, DETAP,  /*GRID  PT  SPACING  IN  ETA  */ 

KDE,  /^THERMAL  COND  X  D  ETA  AVE  */ 

CONVECTN, 

POHMIC , 

PB.PBP,  /^ENERGY  STORED  IN  EXCITED  STAT*/ 

SIGMA) 

FLOAT  DEC (16), 

/*TEMPORARY  DENSITY  EQUATION  VARIABLES  */ 

(GAMMAM , GAMMAP , 

D1B,D21B,D32B,P0,IB,NUE,  /*PARAMETERS  FOR  MSOURCE  */ 

A, 

NES2) 

FLOAT  DEC (16), 

(J,IJ)  FIXED  BIN(31) ; 

ON  FINISH  PUT  SKIP(5)  DATA; 

/*SET  THERMAL  &  ELECTRICAL  CONDUCTIVITIES  AT  0+  (E)  &,  1-  (C) .  */ 
IF  TAU ( 1 ) <0 . 1  THEN  DO;  TAU(1)=0.1;  EFIX-’l’B;  END; 

IF  TAU(N)<0. 1  THEN  DO;  TAU(N)=0.1;  CFIX^’l’B;  END; 

IF  RE=0 . 5  THEN  MUEA=TAUN; 

ELSE  IF  RE=0 . 0  THEN  MUEA=TAUN/SQRT(TAU) ; 

ELSE  IF  RE=- .5  THEN  MUEA=TAUN/TAU ; 

ELSE  MUEA=TAUN* (  TAU** (RE -0.5)  ); 

K=(  (RE+2)/FYEN  )*MUEA*NEB*TAU ; 

PC=NEB* (TAU+TAUN) ; 

DETA ,DETAP=LOG (K(2)/K(l))  *  DZ/ (K(2)-K(l)) ; 

/*DETERMINE  EMITTER  SHEATH  */ 

CALL  SHEATH (ECHI ,ENE,I*KN/NEB(1) ,TAU(1) ,ENR/NEB(1) , 

TE,  SMR,  EALPHA,  0.8,  ERR); 

IF  ECHI<=lE-5  |  ECHI>=20  THEN  EFIX='l’B; 

/*FIND  EMITTER  (0+)  DERIVATIVES  FROM  B.C.  */ 

ET_ETA= (TAU ( 1 ) - 1 )*ENE*NEB ( 1 ) /KN  -  I*(ECHI-TAU(l)/2) ; 
ET_Z=ET_ETA/K ( 1 ) ; 

EPC_Z= (SQRT(PI/ 8/EALPHA) /LAMDAR/KN )*NEB ( 1 ) /MUI ( 1 ) 

-  I/MUEA(1) ; 

EN_Z*(  EPC_Z -NEB ( 1 ) * (ET  Z+DTAUNDZ)  )/(  TAU(1)+TAUN(1)  ); 

/*SOLVE  COLLECTOR  SHEATH  “  */ 

CALPHA=1/TAU(N) ; 

CNE=0; 

U»1 . 0 ; GU=G (U) ; DELU= .01; 

CON: DO  IJ=1  TO  50; 

DGDU= (G ( U+DELU )-GU)/DELU; 

DELTAU=-GU/DGDU ; 

U=U+DELTAU; 

GU=G(U) ; 
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IF ( ABS (GU) <= .0001)  THEN  GO  TO  FIN; 

END  CON; 

PUT  SKIP  LIST('  COLLECTOR  SHEATH  FAILED  TO  CONVERGE'); 

STOP; 

CCHI=0;  CFIX='l'B; 

GO  TO  FIN2 ; 

G:PROC(UX)  RETURNS (DEC( 16)); 

DCL(UX,GX,SQX)  DEC(16); 

IF  (UX<=0.)  THEN  SQX=0 . ; 

ELSE  SQX=SQRT(UX) ; 

GX=EXP(UX)*(1.+ERF(SQX))*2. 

-NEB(N)*SQRT(TAU(N))/I/KN; 

RETURN (GX); 

END  G; 

FIN:CCHI=U*TAU(N); 

FIN2 : 

/*DETERMINE  DERIVATIVES  AT  COLLECTOR  (1-)  FROM  B.C. 

CT  ETA=-I*(CCHI-TAU(N)/2); 

CT~Z=CT_ETA/K(N) ; 

CPC  Z=(-SQRT(PI/8/CALPHA)/LAMDAR/KN)*NEB(N)/MUI(N) 

-  I/MUEA(N); 

CN_Z=(  CPC_Z-NEB(N)*(CT_Z+DTAUNDZ)  )  /  C  TAU(N)+TAUN(N)  ); 
CDETA=L0G(K(N)/K(N-1))  *  DZ/(K(N)-K(N-1)) ; 

NBCOA=TAU ( 1 )+TAUN ( 1 ) ; 

NBCOB=SQRT(PI /EALPHA/8 . ) /LAMDAR/KN/MUI ( 1 ) -ENE*NEB (1)/K(1)* 
(TAU ( 1 ) - 1 . ) -DTAUNDZ ; 

NBCOC=I*NEB ( 1 ) /K ( 1 )* (ECHI -TAU ( 1 ) /2 . ) - I/MUEA ( 1 ) ; 
NBC1A=TAU(N)+TAUN(N) ; 

NBClB--SQRT(PI/CALPHA/8 . ) /LAMDAR/KN/MUI (N)-CNE*NEB(N)/K(N)* 
(TAU(N)-l . )+DTAUNDZ; 

NBC1C=-I*NEB(N)/K(N)*(CCHI-TAU(N)/2.)~I/MUEA(N); 

TBCOA=l . ; 

TBCOB=ENE*NEB ( 1 ) /KN+I/2 . ; 

TBCOC=-ENE*NEB ( 1 ) /KN-I*ECHI ; 

TBC1A=-1. ; 

TBClB=CNE*NEB(N)/KN-I/2. ; 

TBC1C=-CNE*NEB(N)/KN+I*CCHI ; 

NEB A ( 0 ) = (NBCO A*LAMNEB ) / ( 2 . *DZ ) ; 

NEBB(0)=-LAMNEB*NBC0B; 

NEBC(0)=-NBC0A*LAMNEB/(2.*DZ); 

NEBV(0)-NBC0B* ( 1 . -LAMNEB)*NEB ( 1 )+NBC0C- ( 1 . -LAMNEB)* 
NBC0A*(NEB(2)-NEB(0))/(2.*DZ); 
TAUA(0)=TBC0A*LAMTAU/(2.*DETAP) ; 

TAUB ( 0 ) = - LAMTAU*TBCOB ; 

TAUC ( 0 ) =-TBCOA*LAMTAU/ ( 2 . *DETAP ) ; 

TAUV(0)=TBC0B*(1 . -LAMTAU)*TAU(1)+TBC0C- (1 . -LAMTAU)* 

TBCOA* (TAU (2 ) -TAU (0 ) ) / ( 2 . *DETAP) ; 

NEBA (N+ 1 )* (NBC 1A*LAMNEB ) / ( 2 . *DZ ) ; 

NEBB(N+1)=-LAMNEB*NBC1B; 

NEBC(N+1)=-NBC1A*LAMNEB/(2.*DZ); 

NEBV(N+1)=NBC1B*(1 . -LAMNEB )*NEB (N)+NBC1C- (1 . -LAMNEB)* 
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NBC1A*(NEB(N+1)-NEB(N-1))/(2.*DZ); 

TAUA (N+ 1 ) =TBC 1A*LAMTAU / (2 .*CDETA) ; 

TAUB (N+l )=-LAMTAU*TBClB ; 

TAUC (N+l )®-TBC 1A*LAMTAU/ ( 2 . *CDETA) ; 

TAUV (Nfl )=TBC1B* ( 1 . -LAMTAU)*TAU (N)+TBC 1C - ( 1 . -LAMTAU)* 

TBC 1 A* (TAU (N+ 1 ) -TAU (N-l))/(2. *CDETA) ; 

/^INITIALIZE  GAMMAP  &  QKP  FOR  LOOP.  */ 

MNS=MUI(1)*NEB(1)+MUI(2)*NEB(2); 

MURS=MUI(2)/MUEA(2)  +  MUI(1)/MUEA(1)*(  1  -2*DZ*( 

(0.5 -RE )*ET_Z/TAU ( 1 )  -  0.5*DTAUNDZ/TAUN(1))  ); 
GAMMAP=  0.5*(  (  (MUKD+MUIC  2  ))*(PC(  1  )-PC(0)) 

)/DZ  +I*MURS  ); 

QKP® (TAU ( 1 ) -TAU (0 ) ) /DETA ; 

MUIS=MUI (1)+MUI (2) ; 

DO  J=1  TO  N; 

/*UPDATE  FOR  NEW  J.  */ 

GAMMAM=GAMMAP ; 

QKM=QKP ; 

DETA=DETAP ; 

MNSOLD=MNS ; 

MUISOLD=MUIS ; 

MURSOLD=MURS ; 

IF  J-=N 
THEN  DO; 

DETAP® LOG (K(J+1)/K(J))  *  DZ/(K(J+1)-K(J)) ; 
MNS=MUI(J)*NEB(J)+MUI (J+1)*NEB(J+1) ; 

MUIS=MUI (J)+MUI (J+l) ; 

MURS=MUI (J)/MUEA(J)  +  MUI (J+1)/MUEA(J+1) ; 

END; 

ELSE 

MURS=MURS  +2*DZ* (MUI (N)/MUEA(N) )*( ( . 5-RE)*CT_Z/TAU(N) 

-  0 . 5*DTAUNDZ/TAUN (N) ) ; 

/*FIND  AMBIPOLAR  FLUX  AT  J+l/2.  */ 

GAMMAP=0 . 5* (  (  MUIS*(PC(J+1)-PC(J)) 

)/DZ  +I*MURS  ); 

/*FIND  MASS  SOURCE  AT  J.  */ 

A=CA/TAU( J)**A . 5 ; 

NES2=NNB(J)  *  TAU(J)**1 .5  *  EXP(  CSAHA-EGNDB/TAU(J)  ); 
D21B=F(7)*(1+F(8)/TAU(J)) ; 

D32B=F (2 )*EXP (F ( 3 ) /TAU ( J ) ) ; 

IB=A*NES2*(  1+F(1)/NEB(J)  )/(  1+D21B*(1+D32B/NEB(J))/NEB(J) 

); 

P0=l+(  F(4)/NEB(J)  )*(  1+F(5)/NEB(J)  )/(  1+F(6)/NEB(J)  ); 
NUE=NEB(J)*NEB(J)/NES2; 

MSOURCE(J)®NEB(J)*IB*(  1-P0*NUE  ); 

IF(IDEN=1)  THEN 

MSOURCE ( J )=NEB ( J)*A*NES2 ; 

NEBDOT(J)=RMUR*(GAMMAP-GAMMAM)/DZ  +  MSOURCE(J); 

NA ( J ) ®RMUR*MU I S* (TAU ( J+ 1 )+TAUN ( J+ 1 ) ) / 2 . /DZ**2 ; 

NB ( J ) ®RMUR* (MUI S+MUI SOLD )* (TAU ( J ) +TAUN ( J ) ) / 2 . /DZ**2 ; 
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NC ( J)=RMUR*MUISOLD* (TAU ( J- 1 )+TAUN ( J- 1 ) ) /2 . /DZ**2 ; 
ND(J)=I*(MURS-MURS0LD)*RMUR/DZ/2 . ; 

NS(J)=IB*(1. -PO*NUE) ; 

IF(IDEN=1)  THEN 
NS(J)=A*NES2; 

NEBA(J)=-DT*NA(J)*LAMNEB; 

NEBB (J)=l . +DT*NB (J)*LAMNEB -DT*NS ( J)*LAMNEB ; 

NEBC (J)=-DT*NC (J)*LAMNEB ; 

NEBV(J)=NEB(J)+DT*NA(J)*(1 . -LAMNEB)*NEB(J+1)-DT*NB(J)* 

(1. -LAMNEB)*NEB(J)+DT*NC(J)*(1. -LAMNEB)*NEB(J-1)+ 

DT*ND( J)+DT*NS ( J)*( 1 . -LAMNEB )*NEB ( J) ; 

KDE=K ( J ) * (DETA+DETAP ) / 2 ; 

QKP= (TAU ( J+l ) -TAU (J) ) /DETAP ; 

CONVECTN=- ( 1 . 5 )*I* ( DETA*QKP+DETAP*QKM ) / (2*KDE) ; 
SIGMA=NEB(J)*MUEA(J) ; 

POHMIC=I*(  I/SIGMA  +  TAU(J)*(  NEB (J+l) -NEB (J-l)  ) 

/  (2*DZ*NEB(J))) ; 

PBP=(  F(9)*NNR/NR  )*EXP(  -F(10)/TAU(J)  ); 

PB=(  F(11)*NNR/NR  )*EXP(  -F(12)/TAU(J)  ); 

CV(J)=1.5*NEB(J)  +  NNB ( J)* (F ( 10)*PBP+F ( 12 )*PB*NUE) 

/  (TAU(J)*TAU(J) ) ; 

ESOURCE ( J)=-ELOSSB*MSOURCE ( J) 

-  NNB(J)*PB*(  2*NUE*NEBDOT(J)/NEB(J)  ); 

TAUDOT(J)*(  (QKP-QKM)/KDE  +  CONVECTN  +  POHMIC  +  ESOURCE(J)  ) 
/  CV(J); 

TA(J)=1 . / (DETAP*KDE*CV( J) ) ; 

TB ( J)=( 1 . /DETAP+1 . /DETA)/KDE/CV( J) ; 

TC(J)=1./DETA/KDE/CV(J); 

TS ( J)=(CONVECTN+POHMIC+ESOURCE ( J) ) /CV ( J) ; 

TAUA(J)=-DT*LAMTAU*TA(J) ; 

TAUB ( J)=l . +DT*LAMTAU*TB ( J) ; 

TAUC(J)*-DT*LAMTAU*TC (J) ; 

TAUV(J)=TAU(J)+DT*(1 . -LAMTAU)*TA(J)*TAU(J+1) 

-DT* ( 1 . -LAMTAU)*TB (J)*TAU(J) 

+DT* ( 1 . - LAMTAU ) *TC ( J ) *TAU ( J - 1 ) 

+TS(J)*DT; 


IF  CHKDOT  THEN  PUT  SKIP(2)  DATA(NEB(J) ,TAU(J) ,MSOURCE(J) , 
J,PB ,PBP,A, 

D2 IB , D3 2B , PO , I B , NUE , NES2 , QKP , GAMMAP , DETAP , MURS ) ; 

END; 


NEBC(0)=NEBC(0)-NEBC(1)*NEBA(0)/NEBA(1) 
NEBB (0 )=NEBB (0 ) -NEBB ( 1 )*NEBA (0 ) /NEBA ( 1 ) 
NEBV ( 0 )»NEBV ( 0 ) -NEBV ( 1 ) *NEB A ( 0 ) /NEBA ( 1 ) 
NEBA(0)=NEBB(0); 

NEBB(0)=NEBC(0); 
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NEBA (N+l )=NEBA (N+l ) -NEBA (N)*NEBC (N+l ) /NEBC (N) ; 

NEBB (N+l )=NEBB (N+l ) -NEBB (N)*NEBC (N+l ) /NEBC (N) ; 

NEBV (N+l )=NEBV (N+l ) -NEBV (N)*NEBC (N+l) /NEBC (N) ; 

NEBC (N+l )=NEBB (N+l); 

NEBB (N+l )=NEBA (N+l); 

TAUC  (  0  )  =TAUC  (  0  )  -TAUC  ( 1 )  *TAUA  (  0  )  /TAUA  ( 1 ) ; 

TAUB ( 0 ) =TAUB ( 0 ) -TAUB ( 1 ) *TAUA ( 0 ) /TAUA ( 1 ) ; 

TAUV ( 0 ) =TAUV (0 ) -TAUV ( 1 ) *TAUA ( 0 ) /TAUA ( 1 ) ; 

TAUA(0)=TAUB(0); 

TAUB(0)=TAUC(0); 

TAUA ( N+ 1 ) =TAUA (N+l ) -TAUA (N)*TAUC (N+l ) /TAUC (N) ; 

TAUB (N+l )=TAUB (N+l ) -TAUB (N)*TAUC (N+l) /TAUC (N) ; 

TAUV (N+l )=TAUV (N+l ) -TAUV (N)*TAUC (N+l ) /TAUC (N) ; 

TAUC (N+ 1 ) =TAUB ( N+ 1 ) ; 

TAUB (N+1)=TAUA (N+l) ; 

ALPHA (0)=-NEBA(0) /NEBB (0) ; 

LOOP 1: DO  J=1  TO  N; 

ALPHA(J)=-NEBA(J)/(NEBC(J)*ALPHA(J-1)+NEBB(J)); 

END  LOOP1 ; 

ALPHA (N+l )=0. ; 

BETA (0)=NEBV(0) /NEBB (0) ; 

L00P2 : DO  J=1  TO  (N+l); 

BETA(J)=(NEBV(J)-NEBC(J)*BETA(J-1))/(NEBC(J)*ALPHA(J-1)+NEBB(J)); 
END  L00P2 ; 

NEBU(N+1)=BETA(N+1) ; 

L00P3 : DO  J=N  TO  0  BY(-l); 

NEBU ( J ) =ALPHA ( J ) *NEBU (J+1)+BETA(J); 

END  L00P3 ; 

ALPHA ( 0 ) =-TAUA ( 0 ) /TAUB ( 0 ) ; 

L00P4 : DO  J=1  TO  N; 

ALPHA (J )=-TAUA ( J) / (TAUC (J)*ALPHA ( J- 1 )+TAUB (J) ) ; 

END  L00P4 ; 

ALPHA (N+1)=0. ; 

BETA ( 0 ) =TAUV ( 0 ) /TAUB ( 0 ) ; 

LOOPS: DO  J=1  TO  (N+l); 

BETA (J)=(TAUV(J) -TAUC (J )*BETA( J-l) )/ (TAUC (J)*ALPHA(J-1)+TAUB(J)) ; 
END  LOOPS; 

TAUU(N+1)=BETA(N+1); 

LOOP6 : DO  J=N  TO  0  BY(-l); 

TAUU ( J ) =ALPHA ( J ) +TAUU ( J+ 1 ) +BETA ( J ) ; 

END  L00P6 ; 

NEBD0T= (NEBU -NEB ) /DT ; 

TAUD0T= (TAUU-TAU) /DT ; 

IF  CHKDOT  THEN  PUT  SKIP(3)  DATA; 

END  DOT; 

END  PREDCOR; 


►  . 

* 

< 

_  -  -1 

§ 
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C.2  SHEATH:  Computes  the  Sheath  Model  Solution 

■  •■/! 

.  'i 

• 

;  H 

•  ^ 

SHEATH : PROC (OUTCHI , OUTNE , INCUR , INTAU , INNRNP , INTE , INMEMI , 

OUTALP , INDUM , INERR ) ; 

• 

DCL( JBAR , INCUR ,CE , INTAU ,NRP , INNRNP ,TE , INTE , SMR, INMEMI ,ERR , 

UZERO , 

j  -  ■  « 

• 

,* 

PHIZ, INDUM, INERR, Cl ,PHIM, OUTCHI , NEB, OUTNE, OUTALP, VI) 

DEC(16) ; 

• 

UZERO=2 . 1 ; 

CE=INTAU ; 

JBAR=INCUR*2 . /SQRT(3 . 14159 )*SQRT(2 . /CE) ; 

9 

NRP=INNRNP; 

TE=INTE; 

SMR=INMEMI ;  /*  SQUARE  ROOT  OF  MASS  RATIO  */ 

9 

ERR=INERR; 

CI=1 . ; 

ir 

PHIM=REPHIM (JBAR , CE , Cl ) ; 

NEB=RENEB (JBAR , PHIM , CE , Cl ) ; 

PHIZ=LOG (INNRNP /NEB) ; 

VI— REVI (UZERO , PHI Z , PHIM , CE , Cl ) ; 

IF(VI<=.001)  THEN  VI=.001; 

• 

!• 

OUTCHI=PHIM; 

i 

OUTNE =NEB ; 

OUTALP=l . /CE/ (VI )**2 ; 

* 

r 

RETURN; 

• 

REPHIM: PROC (TJBAR.TCE, TCI)  RETURNS (  DEC(16)); 

i  . 

DCL(TJBAR,TCE ,TCI ,TPHIM,TA,TB)  DEC(16); 

TA=( . 5000*SQRT(2 . ) - . 2900*SQRT( 1 . ) ) / 

(SQRT(2. )-SQRT(l. )) ; 

TB=(.5000-.2900)/(SQFT(l./l-)-SQRT(l./2.))i 

TPHIM=(TB/ (TA-TJBAR)  )**2  *TCE/2.1; 

‘ 

j<" 

RETURN (TPHIM) ; 

END  REPHIM; 

i 

* 

RENEB: PROC (TJBAR, TPHIM, TCE, TCI)  RETURNS (  DEC(16)); 

DCL (TJBAR.TCE ,TCI ,TNEB .TPHIM , SQRTPI , SQRTCE , SQPHCE , 

C1.C2.EXPHCE)  DEC(16) ; 

( 

SQRTPI=SQRT(3 . 1415926) ; 

i 

SQRTCE=SQRT(TCE) ; 

SQPHCE=SQRT (TPHIM/TCE ) ; 

C1=1.+ERF(SQPHCE); 

.V* 

■ 

IF(TPHIM>10 . )  THEN 

c 

C2=0 . ; 

ELSE 

1 _ 

r 

C2=EXP (TPHIM)* (  1.  -  ERF(SQRT(TPHIM))  ); 

r 

r 

EXPHCE*EXP( -TPHIM/TCE) ; 

L 

1 

./  •  *'.  •*,  **••*.  •*. 

•*.«**  •*.  •*  •*.  ■*.  •%**.  **V'N.  •“v*  .  •*.  4  ‘  v*  * .  •*. 
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TNEB= (TJBAR*SQRT(TCE/2 . ) +SQRTCE /C 1*EXPHCE / SQRTP I ) / 

( 1 . / SQRTPI+SQRTCE / SQRTPI*C2/C 1*EXPHCE ) ; 

RETURN CTNEB); 

END  RENEB; 

REVI : PROC (TUZERO , TPHI Z , TPH IM , TCE , TC I )  RETURNS (  DEC(16)); 
DCL (TUZERO , TPHI Z , TPHIM , TCE , TCI , TVI , R )  DEC ( 1 6 ) ; 

R=0 . ; 

•  IF(TPHIZ-TPHIM>0.)  THEN  R=SQRT( (TPHI Z-TPHIM) /TCI) ; 


END 

FILMED 
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